• 362,569

# Posts Tagged ‘harmonic oscillator’

## One atom basis phonons in 2D

Posted by peeterjoot on January 19, 2014

[Click here for a PDF of this post with nicer formatting]

Let’s tackle a problem like the 2D problem of the final exam, but first more generally. Instead of a square lattice consider the lattice with the geometry illustrated in fig. 1.1.

Fig 1.1: Oblique one atom basis

Here, $\mathbf{a}$ and $\mathbf{b}$ are the vector differences between the equilibrium positions separating the masses along the $K_1$ and $K_2$ interaction directions respectively. The equilibrium spacing for the cross coupling harmonic forces are

\begin{aligned}\begin{aligned}\mathbf{r} &= (\mathbf{b} + \mathbf{a})/2 \\ \mathbf{s} &= (\mathbf{b} - \mathbf{a})/2.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.1)

Based on previous calculations, we can write the equations of motion by inspection

\begin{aligned}\begin{aligned}m \dot{d}{\mathbf{u}}_\mathbf{n} = &-K_1 \text{Proj}_{\hat{\mathbf{a}}} \sum_\pm \left( { \mathbf{u}_\mathbf{n} - \mathbf{u}_{\mathbf{n} \pm(1, 0)}} \right)^2 \\ &-K_2 \text{Proj}_{\hat{\mathbf{b}}} \sum_\pm \left( { \mathbf{u}_\mathbf{n} - \mathbf{u}_{\mathbf{n} \pm(0, 1)}} \right)^2 \\ &-K_3 \text{Proj}_{\hat{\mathbf{r}}} \sum_\pm \left( { \mathbf{u}_\mathbf{n} - \mathbf{u}_{\mathbf{n} \pm(1, 1)}} \right)^2 \\ &-K_4 \text{Proj}_{\hat{\mathbf{s}}} \sum_\pm \left( { \mathbf{u}_\mathbf{n} - \mathbf{u}_{\mathbf{n} \pm(1, -1)}} \right)^2.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.2)

Inserting the trial solution

\begin{aligned}\mathbf{u}_\mathbf{n} = \frac{1}{{\sqrt{m}}} \boldsymbol{\epsilon}(\mathbf{q}) e^{i( \mathbf{r}_\mathbf{n} \cdot \mathbf{q} - \omega t) },\end{aligned} \hspace{\stretch{1}}(1.3)

and using the matrix form for the projection operators, we have

\begin{aligned}\begin{aligned}\omega^2 \boldsymbol{\epsilon} &=\frac{K_1}{m} \hat{\mathbf{a}} \hat{\mathbf{a}}^\text{T} \boldsymbol{\epsilon}\sum_\pm\left( { 1 - e^{\pm i \mathbf{a} \cdot \mathbf{q}} } \right) \\ & +\frac{K_2}{m} \hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} \boldsymbol{\epsilon}\sum_\pm\left( { 1 - e^{\pm i \mathbf{b} \cdot \mathbf{q}} } \right) \\ & +\frac{K_3}{m} \hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} \boldsymbol{\epsilon}\sum_\pm\left( { 1 - e^{\pm i (\mathbf{b} + \mathbf{a}) \cdot \mathbf{q}} } \right) \\ & +\frac{K_3}{m} \hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} \boldsymbol{\epsilon}\sum_\pm\left( { 1 - e^{\pm i (\mathbf{b} - \mathbf{a}) \cdot \mathbf{q}} } \right) \\ &=\frac{4 K_1}{m} \hat{\mathbf{a}} \hat{\mathbf{a}}^\text{T} \boldsymbol{\epsilon} \sin^2\left( { \mathbf{a} \cdot \mathbf{q}/2 } \right)+\frac{4 K_2}{m} \hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} \boldsymbol{\epsilon} \sin^2\left( { \mathbf{b} \cdot \mathbf{q}/2 } \right) \\ &+\frac{4 K_3}{m} \hat{\mathbf{r}} \hat{\mathbf{r}}^\text{T} \boldsymbol{\epsilon} \sin^2\left( { (\mathbf{b} + \mathbf{a}) \cdot \mathbf{q}/2 } \right)+\frac{4 K_4}{m} \hat{\mathbf{s}} \hat{\mathbf{s}}^\text{T} \boldsymbol{\epsilon} \sin^2\left( { (\mathbf{b} - \mathbf{a}) \cdot \mathbf{q}/2 } \right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.4)

This fully specifies our eigenvalue problem. Writing

\begin{aligned}\begin{aligned}S_1 &= \sin^2\left( { \mathbf{a} \cdot \mathbf{q}/2 } \right) \\ S_2 &= \sin^2\left( { \mathbf{b} \cdot \mathbf{q}/2 } \right) \\ S_3 &= \sin^2\left( { (\mathbf{b} + \mathbf{a}) \cdot \mathbf{q}/2 } \right) \\ S_4 &= \sin^2\left( { (\mathbf{b} - \mathbf{a}) \cdot \mathbf{q}/2 } \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.5.5)

\begin{aligned}\boxed{A = \frac{4}{m}\left( { K_1 S_1 \hat{\mathbf{a}} \hat{\mathbf{a}}^\text{T} + K_2 S_2 \hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} + K_3 S_3 \hat{\mathbf{r}} \hat{\mathbf{r}}^\text{T} + K_4 S_4 \hat{\mathbf{s}} \hat{\mathbf{s}}^\text{T}} \right),}\end{aligned} \hspace{\stretch{1}}(1.0.5.5)

we wish to solve

\begin{aligned}A \boldsymbol{\epsilon} = \omega^2 \boldsymbol{\epsilon} = \lambda \boldsymbol{\epsilon}.\end{aligned} \hspace{\stretch{1}}(1.0.6)

Neglecting the specifics of the matrix at hand, consider a generic two by two matrix

\begin{aligned}A = \begin{bmatrix}a & b \\ c & d\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.6)

for which the characteristic equation is

\begin{aligned}0 &= \begin{vmatrix}\lambda - a & - b \\ -c & \lambda -d \end{vmatrix} \\ &= (\lambda - a)(\lambda - d) - b c \\ &= \lambda^2 - (a + d) \lambda + a d - b c \\ &= \lambda^2 - (Tr A) \lambda + \left\lvert {A} \right\rvert \\ &= \left( {\lambda - \frac{Tr A}{2}} \right)^2- \left( {\frac{Tr A}{2}} \right)^2 + \left\lvert {A} \right\rvert.\end{aligned} \hspace{\stretch{1}}(1.0.6)

So our angular frequencies are given by

\begin{aligned}\omega^2 = \frac{1}{{2}} \left( { Tr A \pm \sqrt{ \left(Tr A\right)^2 - 4 \left\lvert {A} \right\rvert }} \right).\end{aligned} \hspace{\stretch{1}}(1.0.6)

The square root can be simplified slightly

\begin{aligned}\left( {Tr A} \right)^2 - 4 \left\lvert {A} \right\rvert \\ &= (a + d)^2 -4 (a d - b c) \\ &= a^2 + d^2 + 2 a d - 4 a d + 4 b c \\ &= (a - d)^2 + 4 b c,\end{aligned} \hspace{\stretch{1}}(1.0.6)

so that, finally, the dispersion relation is

\begin{aligned}\boxed{\omega^2 = \frac{1}{{2}} \left( { d + a \pm \sqrt{ (d - a)^2 + 4 b c } } \right),}\end{aligned} \hspace{\stretch{1}}(1.0.6)

Our eigenvectors will be given by

\begin{aligned}0 = (\lambda - a) \boldsymbol{\epsilon}_1 - b\boldsymbol{\epsilon}_2,\end{aligned} \hspace{\stretch{1}}(1.0.6)

or

\begin{aligned}\boldsymbol{\epsilon}_1 \propto \frac{b}{\lambda - a}\boldsymbol{\epsilon}_2.\end{aligned} \hspace{\stretch{1}}(1.0.6)

So, our eigenvectors, the vectoral components of our atomic displacements, are

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}b \\ \omega^2 - a\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.6)

or

\begin{aligned}\boxed{\boldsymbol{\epsilon} \propto\begin{bmatrix}2 b \\ d - a \pm \sqrt{ (d - a)^2 + 4 b c }\end{bmatrix}.}\end{aligned} \hspace{\stretch{1}}(1.0.6)

Square lattice

There is not too much to gain by expanding out the projection operators explicitly in general. However, let’s do this for the specific case of a square lattice (as on the exam problem). In that case, our projection operators are

\begin{aligned}\hat{\mathbf{a}} \hat{\mathbf{a}}^\text{T} = \begin{bmatrix}1 \\ 0\end{bmatrix}\begin{bmatrix}1 & 0\end{bmatrix}=\begin{bmatrix}1 & 0 \\ 0 & 0\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.16a)

\begin{aligned}\hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} = \begin{bmatrix}0\\ 1 \end{bmatrix}\begin{bmatrix}0 &1 \end{bmatrix}=\begin{bmatrix}0 & 0 \\ 0 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.16b)

\begin{aligned}\hat{\mathbf{r}} \hat{\mathbf{r}}^\text{T} = \frac{1}{{2}}\begin{bmatrix}1 \\ 1 \end{bmatrix}\begin{bmatrix}1 &1 \end{bmatrix}=\frac{1}{{2}}\begin{bmatrix}1 & 1 \\ 1 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.16c)

\begin{aligned}\hat{\mathbf{s}} \hat{\mathbf{s}}^\text{T} = \frac{1}{{2}}\begin{bmatrix}-1 \\ 1 \end{bmatrix}\begin{bmatrix}-1 &1 \end{bmatrix}=\frac{1}{{2}}\begin{bmatrix}1 & -1 \\ -1 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.16d)

\begin{aligned}\begin{aligned}S_1 &= \sin^2\left( { \mathbf{a} \cdot \mathbf{q} } \right) \\ S_2 &= \sin^2\left( { \mathbf{b} \cdot \mathbf{q} } \right) \\ S_3 &= \sin^2\left( { (\mathbf{b} + \mathbf{a}) \cdot \mathbf{q} } \right) \\ S_4 &= \sin^2\left( { (\mathbf{b} - \mathbf{a}) \cdot \mathbf{q} } \right),\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.16d)

Our matrix is

\begin{aligned}A = \frac{2}{m}\begin{bmatrix}2 K_1 S_1 + K_3 S_3 + K_4 S_4 & K_3 S_3 - K_4 S_4 \\ K_3 S_3 - K_4 S_4 & 2 K_2 S_2 + K_3 S_3 + K_4 S_4\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.16d)

where, specifically, the squared sines for this geometry are

\begin{aligned}S_1 = \sin^2 \left( { \mathbf{a} \cdot \mathbf{q}/2 } \right) = \sin^2 \left( { a q_x/2} \right)\end{aligned} \hspace{\stretch{1}}(1.0.19a)

\begin{aligned}S_2 = \sin^2 \left( { \mathbf{b} \cdot \mathbf{q}/2 } \right) = \sin^2 \left( { a q_y/2} \right)\end{aligned} \hspace{\stretch{1}}(1.0.19b)

\begin{aligned}S_3 = \sin^2 \left( { (\mathbf{b} + \mathbf{a}) \cdot \mathbf{q}/2 } \right) = \sin^2 \left( { a (q_x + q_y)/2} \right)\end{aligned} \hspace{\stretch{1}}(1.0.19c)

\begin{aligned}S_4 = \sin^2 \left( { (\mathbf{b} - \mathbf{a}) \cdot \mathbf{q}/2 } \right) = \sin^2 \left( { a (q_y - q_x)/2} \right).\end{aligned} \hspace{\stretch{1}}(1.0.19d)

Using eq. 1.0.6, the dispersion relation and eigenvectors are

\begin{aligned}\omega^2 = \frac{2}{m} \left( { \sum_i K_i S_i \pm \sqrt{ (K_2 S_2 - K_1 S_1)^2 + (K_3 S_3 - K_4 S_4)^2 } } \right)\end{aligned} \hspace{\stretch{1}}(1.0.20.20)

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}K_3 S_3 - K_4 S_4 \\ K_2 S_2 - K_1 S_1 \pm \sqrt{ (K_2 S_2 - K_1 S_1)^2 + (K_3 S_3 - K_4 S_4)^2 } \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.20.20)

This calculation is confirmed in oneAtomBasisPhononSquareLatticeEigensystem.nb. Mathematica calculates an alternate form (equivalent to using a zero dot product for the second row), of

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}K_1 S_1 - K_2 S_2 \pm \sqrt{ (K_2 S_2 - K_1 S_1)^2 + (K_3 S_3 - K_4 S_4)^2 } \\ K_3 S_3 - K_4 S_4 \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.20.20)

Either way, we see that $K_3 S_3 - K_4 S_4 = 0$ leads to only horizontal or vertical motion.

With the exam criteria

In the specific case that we had on the exam where $K_1 = K_2$ and $K_3 = K_4$, these are

\begin{aligned}\omega^2 = \frac{2}{m} \left( { K_1 (S_1 + S_2) + K_3(S_3 + S_4) \pm \sqrt{ K_1^2 (S_2 - S_1)^2 + K_3^2 (S_3 - S_4)^2 } } \right)\end{aligned} \hspace{\stretch{1}}(1.0.22.22)

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}K_3 \left( { S_3 - S_4 } \right) \\ K_1 \left( { (S_1 - S_2) \pm \sqrt{ (S_2 - S_1)^2 + \left( \frac{K_3}{K_1} \right)^2 (S_3 - S_4)^2 } } \right)\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.22.22)

For horizontal and vertical motion we need $S_3 = S_4$, or for a $2 \pi \times \text{integer}$ difference in the absolute values of the sine arguments

\begin{aligned}\pm ( a (q_x + q_y) /2 ) = a (q_y - q_y) /2 + 2 \pi n.\end{aligned} \hspace{\stretch{1}}(1.0.22.22)

That is, one of

\begin{aligned}\begin{aligned}q_x &= \frac{2 \pi}{a} n \\ q_y &= \frac{2 \pi}{a} n\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.22.22)

In the first BZ, that is one of $q_x = 0$ or $q_y = 0$.

System in rotated coordinates

On the exam, where we were asked to solve for motion along the cross directions explicitly, there was a strong hint to consider a rotated (by $\pi/4$) coordinate system.

The rotated the lattice basis vectors are $\mathbf{a} = a \mathbf{e}_1, \mathbf{b} = a \mathbf{e}_2$, and the projection matrices. Writing $\hat{\mathbf{r}} = \mathbf{f}_1$ and $\hat{\mathbf{s}} = \mathbf{f}_2$, where $\mathbf{f}_1 = (\mathbf{e}_1 + \mathbf{e}_2)/\sqrt{2}, \mathbf{f}_2 = (\mathbf{e}_2 - \mathbf{e}_1)/\sqrt{2}$, or $\mathbf{e}_1 = (\mathbf{f}_1 - \mathbf{f}_2)/\sqrt{2}, \mathbf{e}_2 = (\mathbf{f}_1 + \mathbf{f}_2)/\sqrt{2}$. In the $\{\mathbf{f}_1, \mathbf{f}_2\}$ basis the projection matrices are

\begin{aligned}\hat{\mathbf{a}} \hat{\mathbf{a}}^\text{T} = \frac{1}{{2}}\begin{bmatrix}1 \\ -1\end{bmatrix}\begin{bmatrix}1 & -1\end{bmatrix}= \frac{1}{{2}} \begin{bmatrix}1 & -1 \\ -1 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.25a)

\begin{aligned}\hat{\mathbf{b}} \hat{\mathbf{b}}^\text{T} = \frac{1}{{2}}\begin{bmatrix}1 \\ 1\end{bmatrix}\begin{bmatrix}1 & 1\end{bmatrix}= \frac{1}{{2}} \begin{bmatrix}1 & 1 \\ 1 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.25b)

\begin{aligned}\hat{\mathbf{r}} \hat{\mathbf{r}}^\text{T} = \begin{bmatrix}1 & 0 \\ 0 & 0\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.25c)

\begin{aligned}\hat{\mathbf{s}} \hat{\mathbf{s}}^\text{T} = \begin{bmatrix}0 & 0 \\ 0 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.25d)

The dot products that show up in the squared sines are

\begin{aligned}\mathbf{a} \cdot \mathbf{q}=a \frac{1}{{\sqrt{2}}} (\mathbf{f}_1 - \mathbf{f}_2) \cdot (\mathbf{f}_1 k_u + \mathbf{f}_2 k_v)=\frac{a}{\sqrt{2}} (k_u - k_v)\end{aligned} \hspace{\stretch{1}}(1.0.26a)

\begin{aligned}\mathbf{b} \cdot \mathbf{q}=a \frac{1}{{\sqrt{2}}} (\mathbf{f}_1 + \mathbf{f}_2) \cdot (\mathbf{f}_1 k_u + \mathbf{f}_2 k_v)=\frac{a}{\sqrt{2}} (k_u + k_v)\end{aligned} \hspace{\stretch{1}}(1.0.26b)

\begin{aligned}(\mathbf{a} + \mathbf{b}) \cdot \mathbf{q} = \sqrt{2} a k_u \end{aligned} \hspace{\stretch{1}}(1.0.26c)

\begin{aligned}(\mathbf{b} - \mathbf{a}) \cdot \mathbf{q} = \sqrt{2} a k_v \end{aligned} \hspace{\stretch{1}}(1.0.26d)

So that in this basis

\begin{aligned}\begin{aligned}S_1 &= \sin^2 \left( { \frac{a}{\sqrt{2}} (k_u - k_v) } \right) \\ S_2 &= \sin^2 \left( { \frac{a}{\sqrt{2}} (k_u + k_v) } \right) \\ S_3 &= \sin^2 \left( { \sqrt{2} a k_u } \right) \\ S_4 &= \sin^2 \left( { \sqrt{2} a k_v } \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.26d)

With the rotated projection operators eq. 1.0.5.5 takes the form

\begin{aligned}A = \frac{2}{m}\begin{bmatrix}K_1 S_1 + K_2 S_2 + 2 K_3 S_3 & K_2 S_2 - K_1 S_1 \\ K_2 S_2 - K_1 S_1 & K_1 S_1 + K_2 S_2 + 2 K_4 S_4\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.26d)

This clearly differs from eq. 1.0.16d, and results in a different expression for the eigenvectors, but the same as eq. 1.0.20.20 for the angular frequencies.

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}K_2 S_2 - K_1 S_1 \\ K_4 S_4 - K_3 S_3 \mp \sqrt{ (K_2 S_2 - K_1 S_1)^2 + (K_3 S_3 - K_4 S_4)^2 }\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.26d)

or, equivalently

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}K_4 S_4 - K_3 S_3 \mp \sqrt{ (K_2 S_2 - K_1 S_1)^2 + (K_3 S_3 - K_4 S_4)^2 } \\ K_1 S_1 - K_2 S_2 \\ \end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.26d)

For the $K_1 = K_2$ and $K_3 = K_4$ case of the exam, this is

\begin{aligned}\boldsymbol{\epsilon} \propto\begin{bmatrix}K_1 (S_2 - S_1 ) \\ K_3 \left( { S_4 - S_3 \mp \sqrt{ \left( \frac{K_1}{K_3} \right)^2 (S_2 - S_1)^2 + (S_3 - S_4)^2 } } \right)\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.26d)

Similar to the horizontal coordinate system, we see that we have motion along the diagonals when

\begin{aligned}\pm \frac{a}{\sqrt{2}} (k_u - k_v) = \frac{a}{\sqrt{2}} (k_u + k_v) + 2 \pi n,\end{aligned} \hspace{\stretch{1}}(1.0.26d)

or one of

\begin{aligned}\begin{aligned}k_u &= \sqrt{2} \frac{\pi}{a} n \\ k_v &= \sqrt{2} \frac{\pi}{a} n\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.26d)

Stability?

The exam asked why the cross coupling is required for stability. Clearly we have more complex interaction. The constant $\omega$ surfaces will also be more complex. However, I still don’t have a good intuition what exactly was sought after for that part of the question.

Numerical computations

A Manipulate allowing for choice of the spring constants and lattice orientation, as shown in fig. 1.2, is available in phy487/oneAtomBasisPhonon.nb. This interface also provides a numerical calculation of the distribution relation as shown in fig. 1.3, and provides an animation of the normal modes for any given selection of $\mathbf{q}$ and $\omega(\mathbf{q})$ (not shown).

Fig 1.2: 2D Single atom basis Manipulate interface

Fig 1.3: Sample distribution relation for 2D single atom basis.

## Two body harmonic oscillator in 3D, and 2D diamond lattice vibrations

Posted by peeterjoot on January 7, 2014

# Abridged harmonic oscillator notes

[This is an abbreviation of more extensive PDF notes associated with the latter part of this post.]

# Motivation and summary of harmonic oscillator background

After having had some trouble on a non-1D harmonic oscillator lattice problem on the exam, I attempted such a problem with enough time available to consider it properly. I found it helpful to consider first just two masses interacting harmonically in 3D, each displaced from an equilibrium position.

The Lagrangian that described this most naturally was found to be

\begin{aligned}\mathcal{L} = \frac{1}{2} m_1 \left( \dot{\mathbf{r}}_1 \right)^2+\frac{1}{2} m_2 \left( \dot{\mathbf{r}}_2 \right)^2- \frac{K}{2} \left( \left\lvert {\mathbf{r}_2 - \mathbf{r}_1} \right\rvert - a \right)^2.\end{aligned} \hspace{\stretch{1}}(2.1)

This was solved in absolute and displacement coordinates, and then I moved on to consider a linear expansion of the harmonic potential about the equilibrium point, a problem closer to the exam problem (albeit still considering only two masses). The equilibrium points were described with vectors $\mathbf{a}_1, \mathbf{a}_2$ as in fig. 2.1, where $\Delta \mathbf{a} = \left\lvert {\Delta \mathbf{a}} \right\rvert (\cos \theta_1, \cos\theta_2, \cos\theta_3)$.

fig 2.1: Direction cosines relative to equilibrium position difference vector

Using such coordinates, and generalizing, it was found that the complete Lagrangian, to second order about the equilibrium positions, is

\begin{aligned}\mathcal{L} = \sum_j \frac{m_i}{2} \dot{u}_{ij}^2 -\frac{K}{2} \sum_{i j} \cos\theta_i \cos\theta_j \left( u_{2 i} - u_{1 i} \right)\left( u_{2 j} - u_{1 j} \right).\end{aligned} \hspace{\stretch{1}}(2.2)

Evaluating the Euler-Lagrange equations, the equations of motion for the displacements were found to be

\begin{aligned}\begin{aligned}m_1 \ddot{\mathbf{u}}_1 &= K \widehat{\Delta \mathbf{a}} \left( \widehat{\Delta \mathbf{a}} \cdot \Delta \mathbf{u} \right) \\ m_2 \ddot{\mathbf{u}}_2 &= -K \widehat{\Delta \mathbf{a}} \left( \widehat{\Delta \mathbf{a}} \cdot \Delta \mathbf{u} \right),\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.3)

or

\begin{aligned}\boxed{\begin{aligned}\mu \Delta \ddot{\mathbf{u}} &= -K \widehat{\Delta \mathbf{a}} \left( \widehat{\Delta \mathbf{a}} \cdot \Delta \mathbf{u} \right) \\ m_1 \ddot{\mathbf{u}}_1 + m_2 \ddot{\mathbf{u}}_2 &= 0.\end{aligned}}\end{aligned} \hspace{\stretch{1}}(2.4)

Observe that on the RHS above we have a projection operator, so we could also write

\begin{aligned}\mu \Delta \ddot{\mathbf{u}} = -K \text{Proj}_{\widehat{\Delta \mathbf{a}}} \Delta \mathbf{u}.\end{aligned} \hspace{\stretch{1}}(2.5)

We see that the equations of motion for the displacements of a system of harmonic oscillators has a rather pleasant expression in terms of projection operators, where we have projections onto the unit vectors between each pair of equilibrium position.

# A number of harmonically coupled masses

Now let’s consider masses at lattice points indexed by a lattice vector $\mathbf{n}$, as illustrated in fig. 2.2.

fig 2.2: Masses harmonically coupled in a lattice

With a coupling constant of $K_{\mathbf{n} \mathbf{m}}$ between lattice points indexed $\mathbf{n}$ and $\mathbf{m}$ (located at $\mathbf{a}_\mathbf{n}$ and $\mathbf{a}_\mathbf{m}$ respectively), and direction cosines for the equilibrium direction vector between those points given by

\begin{aligned}\mathbf{a}_\mathbf{n} - \mathbf{a}_\mathbf{m} = \Delta \mathbf{a}_{\mathbf{n} \mathbf{m}}= \left\lvert {\Delta \mathbf{a}_{\mathbf{n} \mathbf{m}}} \right\rvert (\cos \theta_{\mathbf{n} \mathbf{m} 1},\cos \theta_{\mathbf{n} \mathbf{m} 2},\cos \theta_{\mathbf{n} \mathbf{m} 3}),\end{aligned} \hspace{\stretch{1}}(2.6)

the Lagrangian is

\begin{aligned}\mathcal{L} = \sum_{\mathbf{n}, i} \frac{m_\mathbf{n}}{2} \dot{u}_{\mathbf{n} i}^2-\frac{1}{2} \sum_{\mathbf{n} \ne \mathbf{m}, i, j} \frac{K_{\mathbf{n} \mathbf{m}}}{2} \cos\theta_{\mathbf{n} \mathbf{m} i}\cos\theta_{\mathbf{n} \mathbf{m} j}\left( u_{\mathbf{n} i} - u_{\mathbf{m} i} \right)\left( u_{\mathbf{n} j} - u_{\mathbf{m} j} \right)\end{aligned} \hspace{\stretch{1}}(2.7)

Evaluating the Euler-Lagrange equations for the mass at index $\mathbf{n}$ we have

\begin{aligned}\frac{d}{dt} \frac{\partial {\mathcal{L}}}{\partial {\dot{u}_{\mathbf{n} k}}} =m_\mathbf{n} \ddot{u}_{\mathbf{n} k},\end{aligned} \hspace{\stretch{1}}(2.8)

and

\begin{aligned}\frac{\partial {\mathcal{L}}}{\partial {u_{\mathbf{n} k}}} &= -\sum_{\mathbf{m}, i, j}\frac{K_{\mathbf{n} \mathbf{m}}}{2} \cos\theta_{\mathbf{n} \mathbf{m} i}\cos\theta_{\mathbf{n} \mathbf{m} j}\left(\delta_{i k}\left( u_{\mathbf{n} j} - u_{\mathbf{m} j} \right)+\left( u_{\mathbf{n} i} - u_{\mathbf{m} i} \right)\delta_{j k}\right) \\ &= -\sum_{\mathbf{m}, i}K_{\mathbf{n} \mathbf{m}}\cos\theta_{\mathbf{n} \mathbf{m} k}\cos\theta_{\mathbf{n} \mathbf{m} i}\left( u_{\mathbf{n} i} - u_{\mathbf{m} i} \right) \\ &= -\sum_{\mathbf{m}}K_{\mathbf{n} \mathbf{m}}\cos\theta_{\mathbf{n} \mathbf{m} k}\widehat{\Delta \mathbf{a}} \cdot \Delta \mathbf{u}_{\mathbf{n} \mathbf{m}},\end{aligned} \hspace{\stretch{1}}(2.9)

where $\Delta \mathbf{u}_{\mathbf{n} \mathbf{m}} = \mathbf{u}_\mathbf{n} - \mathbf{u}_\mathbf{m}$. Equating both, we have in vector form

\begin{aligned}m_\mathbf{n} \ddot{\mathbf{u}}_\mathbf{n} = -\sum_{\mathbf{m}}K_{\mathbf{n} \mathbf{m}}\widehat{\Delta \mathbf{a}}\left( \widehat{\Delta \mathbf{a}} \cdot \Delta \mathbf{u}_{\mathbf{n} \mathbf{m}} \right),\end{aligned} \hspace{\stretch{1}}(2.10)

or

\begin{aligned}\boxed{m_\mathbf{n} \ddot{\mathbf{u}}_\mathbf{n} = -\sum_{\mathbf{m}}K_{\mathbf{n} \mathbf{m}}\text{Proj}_{ \widehat{\Delta \mathbf{a}} } \Delta \mathbf{u}_{\mathbf{n} \mathbf{m}},}\end{aligned} \hspace{\stretch{1}}(2.11)

This is an intuitively pleasing result. We have displacement and the direction of the lattice separations in the mix, but not the magnitude of the lattice separation itself.

# Two atom basis, 2D diamond lattice

As a concrete application of the previously calculated equilibrium harmonic oscillator result, let’s consider a two atom basis diamond lattice where the horizontal length is $a$ and vertical height is $b$.

Indexing for the primitive unit cells is illustrated in fig. 2.3.

fig 2.3: Primitive unit cells for rectangular lattice

Let’s write

\begin{aligned}\begin{aligned}\mathbf{r} &= a (\cos\theta, \sin\theta) = a \hat{\mathbf{r}} \\ \mathbf{s} &= a (-\cos\theta, \sin\theta) = a \hat{\mathbf{s}} \\ \mathbf{n} &= (n_1, n_2) \\ \mathbf{r}_\mathbf{n} &= n_1 \mathbf{r} + n_2 \mathbf{s},\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.12)

For mass $m_\alpha, \alpha \in \{1, 2\}$ assume a trial solution of the form

\begin{aligned}\mathbf{u}_{\mathbf{n},\alpha} = \frac{\boldsymbol{\epsilon}_\alpha(\mathbf{q})}{\sqrt{m_\alpha}} e^{i \mathbf{r}_n \cdot \mathbf{q} - \omega t}.\end{aligned} \hspace{\stretch{1}}(2.13)

The equations of motion for the two particles are

\begin{aligned}\begin{aligned}m_1 \ddot{\mathbf{u}}_{\mathbf{n}, 1} &= - K_1 \text{Proj}_{\hat{\mathbf{x}}} \left( \mathbf{u}_{\mathbf{n}, 1} - \mathbf{u}_{\mathbf{n} - (0,1), 2} \right)- K_1 \text{Proj}_{\hat{\mathbf{x}}} \left( \mathbf{u}_{\mathbf{n}, 1} - \mathbf{u}_{\mathbf{n} - (1,0), 2} \right) \\ & \quad- K_2 \text{Proj}_{\hat{\mathbf{y}}} \left( \mathbf{u}_{\mathbf{n}, 1} - \mathbf{u}_{\mathbf{n}, 2} \right)- K_2 \text{Proj}_{\hat{\mathbf{y}}} \left( \mathbf{u}_{\mathbf{n}, 1} - \mathbf{u}_{\mathbf{n} - (1,1), 2} \right) \\ & \quad- K_3 \sum_\pm\text{Proj}_{\hat{\mathbf{r}}} \left( \mathbf{u}_{\mathbf{n}, 1} - \mathbf{u}_{\mathbf{n} \pm (1,0), 1} \right)- K_4 \sum_\pm\text{Proj}_{\hat{\mathbf{s}}} \left( \mathbf{u}_{\mathbf{n}, 1} - \mathbf{u}_{\mathbf{n} \pm (0,1), 1} \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.0.14.14)

\begin{aligned}\begin{aligned}m_2 \ddot{\mathbf{u}}_{\mathbf{n}, 2} &= - K_1 \text{Proj}_{\hat{\mathbf{x}}} \left( \mathbf{u}_{\mathbf{n}, 2} - \mathbf{u}_{\mathbf{n} + (1,0), 1} \right)- K_1 \text{Proj}_{\hat{\mathbf{x}}} \left( \mathbf{u}_{\mathbf{n}, 2} - \mathbf{u}_{\mathbf{n} + (0,1), 1} \right)\\ &\quad- K_2 \text{Proj}_{\hat{\mathbf{y}}} \left( \mathbf{u}_{\mathbf{n}, 2} - \mathbf{u}_{\mathbf{n}, 1} \right)- K_2 \text{Proj}_{\hat{\mathbf{y}}} \left( \mathbf{u}_{\mathbf{n}, 2} - \mathbf{u}_{\mathbf{n} + (1,1), 1} \right)\\ &\quad- K_3 \sum_\pm\text{Proj}_{\hat{\mathbf{r}}} \left( \mathbf{u}_{\mathbf{n}, 2} - \mathbf{u}_{\mathbf{n} \pm (1,0), 2} \right)- K_4 \sum_\pm\text{Proj}_{\hat{\mathbf{s}}} \left( \mathbf{u}_{\mathbf{n}, 2} - \mathbf{u}_{\mathbf{n} \pm (0,1), 2} \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.0.14.14)

Insertion of the trial solution gives

\begin{aligned}\begin{aligned} \omega^2 \sqrt{m_1} \boldsymbol{\epsilon}_1&= K_1 \text{Proj}_{\hat{\mathbf{x}}} \left( \frac{\boldsymbol{\epsilon}_1}{\sqrt{m_1}} - \frac{\boldsymbol{\epsilon}_2}{\sqrt{m_2}} e^{ - i \mathbf{s} \cdot \mathbf{q} } \right)+ K_1 \text{Proj}_{\hat{\mathbf{x}}} \left( \frac{\boldsymbol{\epsilon}_1}{\sqrt{m_1}} - \frac{\boldsymbol{\epsilon}_2}{\sqrt{m_2}} e^{ - i \mathbf{r} \cdot \mathbf{q} } \right) \\ &\quad+ K_2 \text{Proj}_{\hat{\mathbf{y}}} \left( \frac{\boldsymbol{\epsilon}_1}{\sqrt{m_1}} - \frac{\boldsymbol{\epsilon}_2}{\sqrt{m_2}} \right)+ K_2 \text{Proj}_{\hat{\mathbf{y}}} \left( \frac{\boldsymbol{\epsilon}_1}{\sqrt{m_1}} - \frac{\boldsymbol{\epsilon}_2}{\sqrt{m_2}} e^{ - i (\mathbf{r} + \mathbf{s}) \cdot \mathbf{q} } \right) \\ &\quad+ K_3 \left( \text{Proj}_{\hat{\mathbf{r}}} \frac{\boldsymbol{\epsilon}_1}{\sqrt{m_1}} \right)\sum_\pm\left( 1 - e^{ \pm i \mathbf{r} \cdot \mathbf{q} } \right)+ K_4 \left( \text{Proj}_{\hat{\mathbf{s}}} \frac{\boldsymbol{\epsilon}_1}{\sqrt{m_1}} \right)\sum_\pm\left( 1 - e^{ \pm i \mathbf{s} \cdot \mathbf{q} } \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.0.15.15)

\begin{aligned}\begin{aligned}\omega^2 \sqrt{m_2} \boldsymbol{\epsilon}_2&=K_1 \text{Proj}_{\hat{\mathbf{x}}} \left( \frac{\boldsymbol{\epsilon}_2}{\sqrt{m_2}} - \frac{\boldsymbol{\epsilon}_1}{\sqrt{m_1}} e^{ + i \mathbf{r} \cdot \mathbf{q} } \right)+ K_1 \text{Proj}_{\hat{\mathbf{x}}} \left( \frac{\boldsymbol{\epsilon}_2}{\sqrt{m_2}} - \frac{\boldsymbol{\epsilon}_1}{\sqrt{m_1}} e^{ + i \mathbf{s} \cdot \mathbf{q} } \right)\\ &\quad+ K_2 \text{Proj}_{\hat{\mathbf{y}}} \left( \frac{\boldsymbol{\epsilon}_2}{\sqrt{m_2}} - \frac{\boldsymbol{\epsilon}_1}{\sqrt{m_1}} \right)+ K_2 \text{Proj}_{\hat{\mathbf{y}}} \left( \frac{\boldsymbol{\epsilon}_2}{\sqrt{m_2}} - \frac{\boldsymbol{\epsilon}_1}{\sqrt{m_1}} e^{ + i (\mathbf{r} + \mathbf{s}) \cdot \mathbf{q} } \right) \\ &\quad+ K_3 \left( \text{Proj}_{\hat{\mathbf{r}}} \frac{\boldsymbol{\epsilon}_2}{\sqrt{m_2}} \right)\sum_\pm\left( 1 - e^{ \pm i \mathbf{r} \cdot \mathbf{q} } \right)+ K_4 \left( \text{Proj}_{\hat{\mathbf{s}}} \frac{\boldsymbol{\epsilon}_2}{\sqrt{m_2}} \right)\sum_\pm\left( 1 - e^{ \pm i \mathbf{s} \cdot \mathbf{q} } \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.0.15.15)

Regrouping, and using the matrix form $\text{Proj}_{\hat{\mathbf{u}}} = \hat{\mathbf{u}} \hat{\mathbf{u}}^\text{T}$ for the projection operators, this is

\begin{aligned}\left(\omega^2 - \frac{2}{m_1} \left( K_1 \hat{\mathbf{x}} \hat{\mathbf{x}}^T + K_2 \hat{\mathbf{y}} \hat{\mathbf{y}}^T + 2 K_3 \hat{\mathbf{r}} \hat{\mathbf{r}}^T \sin^2 (\mathbf{r} \cdot \mathbf{q}/2) + 2 K_4 \hat{\mathbf{s}} \hat{\mathbf{s}}^T \sin^2 (\mathbf{s} \cdot \mathbf{q}/2) \right)\right)\boldsymbol{\epsilon}_1 = -\left( K_1 \hat{\mathbf{r}} \hat{\mathbf{r}}^\text{T} \left( e^{ - i \mathbf{s} \cdot \mathbf{q} } + e^{ - i \mathbf{r} \cdot \mathbf{q} } \right) + K_2 \hat{\mathbf{s}} \hat{\mathbf{s}}^\text{T} \left( 1 + e^{ - i (\mathbf{r} + \mathbf{s}) \cdot \mathbf{q} } \right) \right)\frac{\boldsymbol{\epsilon}_2}{\sqrt{ m_1 m_2 }}\end{aligned} \hspace{\stretch{1}}(2.0.16.16)

\begin{aligned}\left( \omega^2 - \frac{2}{m_2} \left( K_1 \hat{\mathbf{x}} \hat{\mathbf{x}}^T + K_2 \hat{\mathbf{y}} \hat{\mathbf{y}}^T + 2 K_3 \hat{\mathbf{r}} \hat{\mathbf{r}}^T \sin^2 (\mathbf{r} \cdot \mathbf{q}/2)+ 2 K_4 \hat{\mathbf{s}} \hat{\mathbf{s}}^T \sin^2 (\mathbf{s} \cdot \mathbf{q}/2) \right) \right)\boldsymbol{\epsilon}_2 = -\left( K_1 \hat{\mathbf{r}} \hat{\mathbf{r}}^\text{T} \left( e^{ i \mathbf{s} \cdot \mathbf{q} } + e^{ i \mathbf{r} \cdot \mathbf{q} } \right) + K_2 \hat{\mathbf{s}} \hat{\mathbf{s}}^\text{T} \left( 1 + e^{ i (\mathbf{r} + \mathbf{s}) \cdot \mathbf{q} } \right) \right)\frac{\boldsymbol{\epsilon}_1}{\sqrt{ m_1 m_2 }}\end{aligned} \hspace{\stretch{1}}(2.0.16.16)

As a single matrix equation, this is

\begin{aligned}A = K_1 \hat{\mathbf{x}} \hat{\mathbf{x}}^T + K_2 \hat{\mathbf{y}} \hat{\mathbf{y}}^T + 2 K_3 \hat{\mathbf{r}} \hat{\mathbf{r}}^T \sin^2 (\mathbf{r} \cdot \mathbf{q}/2)+ 2 K_4 \hat{\mathbf{s}} \hat{\mathbf{s}}^T \sin^2 (\mathbf{s} \cdot \mathbf{q}/2)\end{aligned} \hspace{\stretch{1}}(2.0.17.17)

\begin{aligned}B = e^{ i (\mathbf{r} + \mathbf{s}) \cdot \mathbf{q}/2 }\left( { K_1 \hat{\mathbf{r}} \hat{\mathbf{r}}^\text{T} \cos\left( (\mathbf{r} - \mathbf{s}) \cdot \mathbf{q}/2 \right)+ K_2 \hat{\mathbf{s}} \hat{\mathbf{s}}^\text{T} \cos\left( (\mathbf{r} + \mathbf{s}) \cdot \mathbf{q}/2 \right)} \right)\end{aligned} \hspace{\stretch{1}}(2.0.17.17)

\begin{aligned}0 =\begin{bmatrix}\omega^2 - \frac{2 A}{m_1} & \frac{B^{*}}{\sqrt{m_1 m_2}} \\ \frac{B}{\sqrt{m_1 m_2}} & \omega^2 - \frac{2 A}{m_2} \end{bmatrix}\begin{bmatrix}\boldsymbol{\epsilon}_1 \\ \boldsymbol{\epsilon}_2\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.0.17.17)

Observe that this is an eigenvalue problem $E \mathbf{e} = \omega^2 \mathbf{e}$ for matrix

\begin{aligned}E = \begin{bmatrix}\frac{2 A}{m_1} & -\frac{B^{*}}{\sqrt{m_1 m_2}} \\ -\frac{B}{\sqrt{m_1 m_2}} & \frac{2 A}{m_2} \end{bmatrix},\end{aligned} \hspace{\stretch{1}}(2.0.17.17)

and eigenvalues $\omega^2$.

To be explicit lets put the $A$ and $B$ functions in explicit matrix form. The orthogonal projectors have a simple form

\begin{aligned}\text{Proj}_{\hat{\mathbf{x}}} = \hat{\mathbf{x}} \hat{\mathbf{x}}^\text{T}= \begin{bmatrix}1 \\ 0\end{bmatrix}\begin{bmatrix}1 & 0\end{bmatrix}=\begin{bmatrix}1 & 0 \\ 0 & 0\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.0.19a)

\begin{aligned}\text{Proj}_{\hat{\mathbf{y}}} = \hat{\mathbf{y}} \hat{\mathbf{y}}^\text{T}= \begin{bmatrix}0 \\ 1\end{bmatrix}\begin{bmatrix}0 & 1\end{bmatrix}=\begin{bmatrix}0 & 0 \\ 0 & 1\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.0.19b)

For the $\hat{\mathbf{r}}$ and $\hat{\mathbf{s}}$ projection operators, we can use half angle formulations

\begin{aligned}\text{Proj}_{\hat{\mathbf{r}}} = \hat{\mathbf{r}} \hat{\mathbf{r}}^\text{T}= \begin{bmatrix}\cos\theta \\ \sin\theta\end{bmatrix}\begin{bmatrix}\cos\theta & \sin\theta\end{bmatrix}=\begin{bmatrix}\cos^2\theta & \cos\theta \sin\theta \\ \cos\theta \sin\theta & \sin^2 \theta\end{bmatrix}=\frac{1}{2}\begin{bmatrix}1 + \cos \left( 2 \theta \right) & \sin \left( 2 \theta \right) \\ \sin \left( 2 \theta \right) & 1 - \cos \left( 2 \theta \right)\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.0.20.20)

\begin{aligned}\text{Proj}_{\hat{\mathbf{s}}} = \hat{\mathbf{s}} \hat{\mathbf{s}}^\text{T}= \begin{bmatrix}-\cos\theta \\ \sin\theta\end{bmatrix}\begin{bmatrix}-\cos\theta & \sin\theta\end{bmatrix}=\begin{bmatrix}\cos^2\theta & -\cos\theta \sin\theta \\ -\cos\theta \sin\theta & \sin^2 \theta\end{bmatrix}=\frac{1}{2}\begin{bmatrix}1 + \cos \left( 2 \theta \right) & -\sin \left( 2 \theta \right) \\ -\sin \left( 2 \theta \right) & 1 - \cos \left( 2 \theta \right)\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.0.20.20)

After some manipulation, and the following helper functions

\begin{aligned}\begin{aligned}\alpha_\pm &= K_3 \sin^2 (\mathbf{r} \cdot \mathbf{q}/2) \pm K_4 \sin^2 (\mathbf{s} \cdot \mathbf{q}/2) \\ \beta_\pm &= K_1 \cos\left( (\mathbf{r} - \mathbf{s}) \cdot \mathbf{q}/2 \right) \pm K_2 \cos\left( (\mathbf{r} + \mathbf{s}) \cdot \mathbf{q}/2 \right),\end{aligned}\end{aligned} \hspace{\stretch{1}}(2.0.20.20)

the block matrices of eq. 2.0.17.17 take the form

\begin{aligned}A = \begin{bmatrix}K_1 + \alpha_+ (1 + \cos\left( 2 \theta \right)) & \alpha_- \sin\left( 2 \theta \right) \\ \alpha_- \sin\left( 2 \theta \right) & K_2 + \alpha_+ (1 - \cos\left( 2 \theta \right))\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.0.22.22)

\begin{aligned}B = e^{ i (\mathbf{r} + \mathbf{s}) \cdot \mathbf{q}/2 }\begin{bmatrix} \beta_+ (1 + \cos \left( 2 \theta \right)) & \beta_- \sin \left( 2 \theta \right) \\ \beta_- \sin \left( 2 \theta \right) & \beta_+( 1 -\cos \left( 2 \theta \right))\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(2.0.22.22)

A final bit of simplification for $B$ possible, noting that $\mathbf{r} + \mathbf{s} = 2 a (0, \sin\theta )$, and $\mathbf{r} - \mathbf{s} = 2 a(\cos\theta, 0)$, so

\begin{aligned}\beta_\pm = K_1 \cos\left( a \cos\theta q_x \right) \pm K_2 \cos\left( a \sin\theta q_y \right),\end{aligned} \hspace{\stretch{1}}(2.0.22.22)

and

\begin{aligned}B = e^{ i a \sin\theta q_y }\begin{bmatrix} \beta_+ (1 + \cos \left( 2 \theta \right)) & \beta_- \sin \left( 2 \theta \right) \\ \beta_- \sin \left( 2 \theta \right) & \beta_+( 1 -\cos \left( 2 \theta \right))\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(2.0.22.22)

It isn’t particularly illuminating to expand out the determinant for such a system, even though it can be done symbolically without too much programming. However, what is easy after formulating the matrix for this system, is actually solving it. This is done, and animated, in twoAtomBasisRectangularLatticeDispersionRelation.cdf

## Exponential solutions to second order linear system

Posted by peeterjoot on October 20, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

# Motivation

We’re discussing specific forms to systems of coupled linear differential equations, such as a loop of “spring” connected masses (i.e. atoms interacting with harmonic oscillator potentials) as sketched in fig. 1.1.

Fig 1.1: Three springs loop

Instead of assuming a solution, let’s see how far we can get attacking this problem systematically.

# Matrix methods

Suppose that we have a set of $N$ masses constrained to a circle interacting with harmonic potentials. The Lagrangian for such a system (using modulo $N$ indexing) is

\begin{aligned}\mathcal{L} = \frac{1}{{2}} \sum_{k = 0}^{2} m_k \dot{u}_k^2 - \frac{1}{{2}} \sum_{k = 0}^2 \kappa_k \left( u_{k+1} - u_k \right)^2.\end{aligned} \hspace{\stretch{1}}(1.1)

The force equations follow directly from the Euler-Lagrange equations

\begin{aligned}0 = \frac{d{{}}}{dt} \frac{\partial {\mathcal{L}}}{\partial {\dot{u}_{n, \alpha}}}- \frac{\partial {\mathcal{L}}}{\partial {u_{n, \alpha}}}.\end{aligned} \hspace{\stretch{1}}(1.2)

For the simple three particle system depicted above, this is

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m_0 \dot{u}_0^2 +\frac{1}{{2}} m_1 \dot{u}_1^2 +\frac{1}{{2}} m_2 \dot{u}_2^2 - \frac{1}{{2}} \kappa_0 \left( u_1 - u_0 \right)^2- \frac{1}{{2}} \kappa_1 \left( u_2 - u_1 \right)^2- \frac{1}{{2}} \kappa_2 \left( u_0 - u_2 \right)^2,\end{aligned} \hspace{\stretch{1}}(1.3)

with equations of motion

\begin{aligned}\begin{aligned}0 &= m_0 \dot{d}{u}_0 + \kappa_0 \left( u_0 - u_1 \right) + \kappa_2 \left( u_0 - u_2 \right) \\ 0 &= m_1 \dot{d}{u}_1 + \kappa_1 \left( u_1 - u_2 \right) + \kappa_0 \left( u_1 - u_0 \right) \\ 0 &= m_2 \dot{d}{u}_2 + \kappa_2 \left( u_2 - u_0 \right) + \kappa_1 \left( u_2 - u_1 \right),\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.4)

Let’s partially non-dimensionalize this. First introduce average mass $\bar{{m}}$ and spring constants $\bar{{\kappa}}$, and rearrange slightly

\begin{aligned}\begin{aligned}\frac{\bar{{m}}}{\bar{{k}}} \dot{d}{u}_0 &= -\frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_0} \left( u_0 - u_1 \right) - \frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_0} \left( u_0 - u_2 \right) \\ \frac{\bar{{m}}}{\bar{{k}}} \dot{d}{u}_1 &= -\frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_1} \left( u_1 - u_2 \right) - \frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_1} \left( u_1 - u_0 \right) \\ \frac{\bar{{m}}}{\bar{{k}}} \dot{d}{u}_2 &= -\frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_2} \left( u_2 - u_0 \right) - \frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_2} \left( u_2 - u_1 \right).\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.5)

With

\begin{aligned}\tau = \sqrt{\frac{\bar{{k}}}{\bar{{m}}}} t = \Omega t\end{aligned} \hspace{\stretch{1}}(1.0.6.6)

\begin{aligned}\mathbf{u} = \begin{bmatrix}u_0 \\ u_1 \\ u_2\end{bmatrix}\end{aligned} \hspace{\stretch{1}}(1.0.6.6)

\begin{aligned}B = \begin{bmatrix}-\frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_0} - \frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_0} &\frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_0} &\frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_0} \\ \frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_1} &-\frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_1} - \frac{\kappa_0 \bar{{m}}}{\bar{{k}} m_1} &\frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_1} \\ \frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_2} & \frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_2} &-\frac{\kappa_2 \bar{{m}}}{\bar{{k}} m_2} - \frac{\kappa_1 \bar{{m}}}{\bar{{k}} m_2} \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.6.6)

Our system takes the form

\begin{aligned}\frac{d^2 \mathbf{u}}{d\tau^2} = B \mathbf{u}.\end{aligned} \hspace{\stretch{1}}(1.0.6.6)

We can at least theoretically solve this in a simple fashion if we first convert it to a first order system. We can do that by augmenting our vector of displacements with their first derivatives

\begin{aligned}\mathbf{w} =\begin{bmatrix}\mathbf{u} \\ \frac{d \mathbf{u}}{d\tau} \end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.6.6)

So that

\begin{aligned}\frac{d \mathbf{w}}{d\tau} =\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} \mathbf{w}= A \mathbf{w}.\end{aligned} \hspace{\stretch{1}}(1.0.9)

Now the solution is conceptually trivial

\begin{aligned}\mathbf{w} = e^{A \tau} \mathbf{w}_\circ.\end{aligned} \hspace{\stretch{1}}(1.0.9)

We are however, faced with the task of exponentiating the matrix $A$. All the powers of this matrix $A$ will be required, but they turn out to be easy to calculate

\begin{aligned}{\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} }^2=\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} \begin{bmatrix}0 & I \\ B & 0\end{bmatrix} =\begin{bmatrix}B & 0 \\ 0 & B\end{bmatrix} \end{aligned} \hspace{\stretch{1}}(1.0.11a)

\begin{aligned}{\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} }^3=\begin{bmatrix}B & 0 \\ 0 & B\end{bmatrix} \begin{bmatrix}0 & I \\ B & 0\end{bmatrix} =\begin{bmatrix}0 & B \\ B^2 & 0\end{bmatrix} \end{aligned} \hspace{\stretch{1}}(1.0.11b)

\begin{aligned}{\begin{bmatrix}0 & I \\ B & 0\end{bmatrix} }^4=\begin{bmatrix}0 & B \\ B^2 & 0\end{bmatrix} \begin{bmatrix}0 & I \\ B & 0\end{bmatrix} =\begin{bmatrix}B^2 & 0 \\ 0 & B^2\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.11c)

allowing us to write out the matrix exponential

\begin{aligned}e^{A \tau} = \sum_{k = 0}^\infty \frac{\tau^{2k}}{(2k)!} \begin{bmatrix}B^k & 0 \\ 0 & B^k\end{bmatrix}+\sum_{k = 0}^\infty \frac{\tau^{2k + 1}}{(2k + 1)!} \begin{bmatrix}0 & B^k \\ B^{k+1} & 0\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.11c)

Case I: No zero eigenvalues

Provided that $B$ has no zero eigenvalues, we could factor this as

\begin{aligned}\begin{bmatrix}0 & B^k \\ B^{k+1} & 0\end{bmatrix}=\begin{bmatrix}0 & B^{-1/2} \\ B^{1/2} & 0\end{bmatrix}\begin{bmatrix}B^{k + 1/2} & 0 \\ 0 & B^{k + 1/2}\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.11c)

This initially leads us to believe the following, but we’ll find out that the three springs interaction matrix $B$ does have a zero eigenvalue, and we’ll have to be more careful. If there were any such interaction matrices that did not have such a zero we could simply write

\begin{aligned}e^{A \tau} = \sum_{k = 0}^\infty \frac{\tau^{2k}}{(2k)!} \begin{bmatrix}\sqrt{B}^{2k} & 0 \\ 0 & \sqrt{B}^{2k}\end{bmatrix}+\begin{bmatrix}0 & B^{-1/2} \\ B^{1/2} & 0\end{bmatrix}\sum_{k = 0}^\infty \frac{\tau^{2k + 1}}{(2k + 1)!} \begin{bmatrix}\sqrt{B}^{2 k + 1} & 0 \\ 0 & \sqrt{B}^{2 k+1} \end{bmatrix}=\cosh\begin{bmatrix}\sqrt{B} \tau & 0 \\ 0 & \sqrt{B} \tau\end{bmatrix}+ \begin{bmatrix}0 & 1/\sqrt{B} \tau \\ \sqrt{B} \tau & 0\end{bmatrix}\sinh\begin{bmatrix}\sqrt{B} \tau & 0 \\ 0 & \sqrt{B} \tau\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.11c)

This is

\begin{aligned}e^{A \tau}=\begin{bmatrix}\cosh \sqrt{B} \tau & (1/\sqrt{B}) \sinh \sqrt{B} \tau \\ \sqrt{B} \sinh \sqrt{B} \tau & \cosh \sqrt{B} \tau \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.11c)

The solution, written out is

\begin{aligned}\begin{bmatrix}\mathbf{u} \\ \mathbf{u}'\end{bmatrix}=\begin{bmatrix}\cosh \sqrt{B} \tau & (1/\sqrt{B}) \sinh \sqrt{B} \tau \\ \sqrt{B} \sinh \sqrt{B} \tau & \cosh \sqrt{B} \tau \end{bmatrix}\begin{bmatrix}\mathbf{u}_\circ \\ \mathbf{u}_\circ'\end{bmatrix},\end{aligned} \hspace{\stretch{1}}(1.0.11c)

so that

\begin{aligned}\boxed{\mathbf{u} = \cosh \sqrt{B} \tau \mathbf{u}_\circ + \frac{1}{{\sqrt{B}}} \sinh \sqrt{B} \tau \mathbf{u}_\circ'.}\end{aligned} \hspace{\stretch{1}}(1.0.11c)

As a check differentiation twice shows that this is in fact the general solution, since we have

\begin{aligned}\mathbf{u}' = \sqrt{B} \sinh \sqrt{B} \tau \mathbf{u}_\circ + \cosh \sqrt{B} \tau \mathbf{u}_\circ',\end{aligned} \hspace{\stretch{1}}(1.0.11c)

and

\begin{aligned}\mathbf{u}'' = B \cosh \sqrt{B} \tau \mathbf{u}_\circ + \sqrt{B} \sinh \sqrt{B} \tau \mathbf{u}_\circ'= B \left( \cosh \sqrt{B} \tau \mathbf{u}_\circ + \frac{1}{{\sqrt{B}}} \sinh \sqrt{B} \tau \mathbf{u}_\circ' \right)= B \mathbf{u}.\end{aligned} \hspace{\stretch{1}}(1.0.11c)

Observe that this solution is a general solution to second order constant coefficient linear systems of the form we have in eq. 1.5. However, to make it meaningful we do have the additional computational task of performing an eigensystem decomposition of the matrix $B$. We expect negative eigenvalues that will give us oscillatory solutions (ie: the matrix square roots will have imaginary eigenvalues).

## Example: An example diagonalization to try things out

{example:threeSpringLoop:1}{

Let’s do that diagonalization for the simplest of the three springs system as an example, with $\kappa_j = \bar{{k}}$ and $m_j = \bar{{m}}$, so that we have

\begin{aligned}B = \begin{bmatrix}-2 & 1 & 1 \\ 1 & -2 & 1 \\ 1 & 1 & -2 \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.20)

A orthonormal eigensystem for $B$ is

\begin{aligned}\left\{\mathbf{e}_{-3, 1}, \mathbf{e}_{-3, 2}, \mathbf{e}_{0, 1} \right\}=\left\{\frac{1}{{\sqrt{6}}}\begin{bmatrix} -1 \\ -1 \\ 2 \end{bmatrix},\frac{1}{{\sqrt{2}}}\begin{bmatrix} -1 \\ 1 \\ 0 \end{bmatrix},\frac{1}{{\sqrt{3}}}\begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}\right\}.\end{aligned} \hspace{\stretch{1}}(1.21)

With

\begin{aligned}\begin{aligned}U &= \frac{1}{{\sqrt{6}}}\begin{bmatrix} -\sqrt{3} & -1 & \sqrt{2} \\ 0 & 2 & \sqrt{2} \\ \sqrt{3} & -1 & \sqrt{2} \end{bmatrix} \\ D &= \begin{bmatrix}-3 & 0 & 0 \\ 0 & -3 & 0 \\ 0 & 0 & 0 \\ \end{bmatrix}\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.22)

We have

\begin{aligned}B = U D U^\text{T},\end{aligned} \hspace{\stretch{1}}(1.23)

We also find that $B$ and its root are intimately related in a surprising way

\begin{aligned}\sqrt{B} = \sqrt{3} iU \begin{bmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \\ \end{bmatrix}U^\text{T}=\frac{1}{{\sqrt{3} i}} B.\end{aligned} \hspace{\stretch{1}}(1.24)

We also see, unfortunately that $B$ has a zero eigenvalue, so we can’t compute $1/\sqrt{B}$. We’ll have to back and up and start again differently.
}

Case II: allowing for zero eigenvalues

Now that we realize we have to deal with zero eigenvalues, a different approach suggests itself. Instead of reducing our system using a Hamiltonian transformation to a first order system, let’s utilize that diagonalization directly. Our system is

\begin{aligned}\mathbf{u}'' = B \mathbf{u} = U D U^{-1} \mathbf{u},\end{aligned} \hspace{\stretch{1}}(1.25)

where $D = [ \lambda_i \delta_{ij} ]$ and

\begin{aligned}\left( U^{-1} \mathbf{u} \right)'' = D \left( U^{-1} \mathbf{u} \right).\end{aligned} \hspace{\stretch{1}}(1.26)

Let

\begin{aligned}\mathbf{z} = U^{-1} \mathbf{u},\end{aligned} \hspace{\stretch{1}}(1.27)

so that our system is just

\begin{aligned}\mathbf{z}'' = D \mathbf{z},\end{aligned} \hspace{\stretch{1}}(1.28)

or

\begin{aligned}z_i'' = \lambda_i z_i.\end{aligned} \hspace{\stretch{1}}(1.29)

This is $N$ equations, each decoupled and solvable by inspection. Suppose we group the eigenvalues into sets $\{ \lambda_n < 0, \lambda_p > 0, \lambda_z = 0 \}$. Our solution is then

\begin{aligned}\mathbf{z} = \sum_{ \lambda_n < 0, \lambda_p > 0, \lambda_z = 0}\left( a_n \cos \sqrt{-\lambda_n} \tau + b_n \sin \sqrt{-\lambda_n} \tau \right)\mathbf{e}_n+\left( a_p \cosh \sqrt{\lambda_p} \tau + b_p \sinh \sqrt{\lambda_p} \tau \right)\mathbf{e}_p+\left( a_z + b_z \tau \right) \mathbf{e}_z.\end{aligned} \hspace{\stretch{1}}(1.30)

Transforming back to lattice coordinates using $\mathbf{u} = U \mathbf{z}$, we have

\begin{aligned}\mathbf{u} = \sum_{ \lambda_n < 0, \lambda_p > 0, \lambda_z = 0}\left( a_n \cos \sqrt{-\lambda_n} \tau + b_n \sin \sqrt{-\lambda_n} \tau \right)U \mathbf{e}_n+\left( a_p \cosh \sqrt{\lambda_p} \tau + b_p \sinh \sqrt{\lambda_p} \tau \right)U \mathbf{e}_p+\left( a_z + b_z \tau \right) U \mathbf{e}_z.\end{aligned} \hspace{\stretch{1}}(1.31)

We see that the zero eigenvalues integration terms have no contribution to the lattice coordinates, since $U \mathbf{e}_z = \lambda_z \mathbf{e}_z = 0$, for all $\lambda_z = 0$.

If $U = [ \mathbf{e}_i ]$ are a set of not necessarily orthonormal eigenvectors for $B$, then the vectors $\mathbf{f}_i$, where $\mathbf{e}_i \cdot \mathbf{f}_j = \delta_{ij}$ are the reciprocal frame vectors. These can be extracted from $U^{-1} = [ \mathbf{f}_i ]^\text{T}$ (i.e., the rows of $U^{-1}$). Taking dot products between $\mathbf{f}_i$ with $\mathbf{u}(0) = \mathbf{u}_\circ$ and $\mathbf{u}'(0) = \mathbf{u}'_\circ$, provides us with the unknown coefficients $a_n, b_n$

\begin{aligned}\mathbf{u}(\tau)= \sum_{ \lambda_n < 0, \lambda_p > 0 }\left( (\mathbf{u}_\circ \cdot \mathbf{f}_n) \cos \sqrt{-\lambda_n} \tau + \frac{\mathbf{u}_\circ' \cdot \mathbf{f}_n}{\sqrt{-\lambda_n} } \sin \sqrt{-\lambda_n} \tau \right)\mathbf{e}_n+\left( (\mathbf{u}_\circ \cdot \mathbf{f}_p) \cosh \sqrt{\lambda_p} \tau + \frac{\mathbf{u}_\circ' \cdot \mathbf{f}_p}{\sqrt{-\lambda_p} } \sinh \sqrt{\lambda_p} \tau \right)\mathbf{e}_p.\end{aligned} \hspace{\stretch{1}}(1.32)

Supposing that we constrain ourself to looking at just the oscillatory solutions (i.e. the lattice does not shake itself to pieces), then we have

\begin{aligned}\boxed{\mathbf{u}(\tau)= \sum_{ \lambda_n < 0 }\left( \sum_j \mathbf{e}_{n,j} \mathbf{f}_{n,j}^\text{T} \right)\left( \mathbf{u}_\circ \cos \sqrt{-\lambda_n} \tau + \frac{\mathbf{u}_\circ' }{\sqrt{-\lambda_n} } \sin \sqrt{-\lambda_n} \tau \right).}\end{aligned} \hspace{\stretch{1}}(1.33)

Eigenvectors for eigenvalues that were degenerate have been explicitly enumerated here, something previously implied. Observe that the dot products of the form $(\mathbf{a} \cdot \mathbf{f}_i) \mathbf{e}_i$ have been put into projector operator form to group terms more nicely. The solution can be thought of as a weighted projector operator working as a time evolution operator from the initial state.

## Example: Our example interaction revisited

{example:threeSpringLoop:2}{

Recall that we had an orthonormal basis for the $\lambda = -3$ eigensubspace for the interaction example of eq. 1.20 again, so $\mathbf{e}_{-3,i} = \mathbf{f}_{-3, i}$. We can sum $\mathbf{e}_{-3,1} \mathbf{e}_{-3, 1}^\text{T} + \mathbf{e}_{-3,2} \mathbf{e}_{-3, 2}^\text{T}$ to find

\begin{aligned}\mathbf{u}(\tau)= \frac{1}{{3}}\begin{bmatrix}2 & -1 & -1 \\ -1 & 2 & -1 \\ -1 & -1 & 2\end{bmatrix}\left( \mathbf{u}_\circ \cos \sqrt{3} \tau + \frac{ \mathbf{u}_\circ' }{\sqrt{3} } \sin \sqrt{3} \tau \right).\end{aligned} \hspace{\stretch{1}}(1.34)

The leading matrix is an orthonormal projector of the initial conditions onto the eigen subspace for $\lambda_n = -3$. Observe that this is proportional to $B$ itself, scaled by the square of the non-zero eigenvalue of $\sqrt{B}$. From this we can confirm by inspection that this is a solution to $\mathbf{u}'' = B \mathbf{u}$, as desired.
}

# Fourier transform methods

Let’s now try another item from our usual toolbox on these sorts of second order systems, the Fourier transform. For a one variable function of time let’s write the transform pair as

\begin{aligned}x(t) = \int_{-\infty}^{\infty} \tilde{x}(\omega) e^{-i \omega t} d\omega\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

\begin{aligned}\tilde{x}(\omega) = \frac{1}{{2\pi}} \int_{-\infty}^{\infty} x(t) e^{i \omega t} dt\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

One mass harmonic oscillator

The simplest second order system is that of the harmonic oscillator

\begin{aligned}0 = \dot{d}{x}(t) + \omega_\circ^2 x.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Application of the transform gives

\begin{aligned}0 = \left( \frac{d^2}{dt^2} + \omega_\circ^2 \right)\int_{-\infty}^{\infty} \tilde{x}(\omega) e^{-i \omega t} d\omega= \int_{-\infty}^{\infty} \left( -\omega^2 + \omega_\circ^2 \right)\tilde{x}(\omega) e^{-i \omega t} d\omega.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

We clearly have a constraint that is a function of frequency, but one that has to hold for all time. Let’s transform this constraint to the frequency domain to consider that constraint independent of time.

\begin{aligned}0 =\frac{1}{{2 \pi}}\int_{-\infty}^{\infty} dt e^{ i \omega t}\int_{-\infty}^{\infty} \left( -{\omega'}^2 + \omega_\circ^2 \right)\tilde{x}(\omega') e^{-i \omega' t} d\omega'=\int_{-\infty}^{\infty} d\omega'\left( -{\omega'}^2 + \omega_\circ^2 \right)\tilde{x}(\omega') \frac{1}{{2 \pi}}\int_{-\infty}^{\infty} dt e^{ i (\omega -\omega') t}=\int_{-\infty}^{\infty} d\omega'\left( -{\omega'}^2 + \omega_\circ^2 \right)\tilde{x}(\omega') \delta( \omega - \omega' )=\left( -{\omega}^2 + \omega_\circ^2 \right)\tilde{x}(\omega).\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

How do we make sense of this?

Since $\omega$ is an integration variable, we can’t just mandate that it equals the constant driving frequency $\pm \omega_\circ$. It’s clear that we require a constraint on the transform $\tilde{x}(\omega)$ as well. As a trial solution, imagine that

\begin{aligned}\tilde{x}(\omega) = \left\{\begin{array}{l l}\tilde{x}_\circ & \quad \mbox{iflatex \left\lvert {\omega – \pm \omega_\circ} \right\rvert < \omega_{\text{cutoff}}} \\ 0 & \quad \mbox{otherwise}\end{array}\right.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

This gives us

\begin{aligned}0 = \tilde{x}_\circ\int_{\pm \omega_\circ -\omega_{\text{cutoff}}}^{\pm \omega_\circ +\omega_{\text{cutoff}}}\left( \omega^2 - \omega_\circ^2 \right)e^{-i \omega t} d\omega.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Now it is clear that we can satisfy our constraint only if the interval $[\pm \omega_\circ -\omega_{\text{cutoff}}, \pm \omega_\circ + \omega_{\text{cutoff}}]$ is made infinitesimal. Specifically, we require both a $\omega^2 = \omega_\circ^2$ constraint and that the transform $\tilde{x}(\omega)$ have a delta function nature. That is

\begin{aligned}\tilde{x}(\omega) = A \delta(\omega - \omega_\circ)+ B \delta(\omega + \omega_\circ).\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Substitution back into our transform gives

\begin{aligned}x(t) = A e^{-i \omega_\circ t}+ B e^{i \omega_\circ t}.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

We can verify quickly that this satisfies our harmonic equation $\dot{d}{x} = -\omega_\circ x$.

Two mass harmonic oscillator

Having applied the transform technique to the very simplest second order system, we can now consider the next more complex system, that of two harmonically interacting masses (i.e. two masses on a frictionless spring).

F1

Our system is described by

\begin{aligned}\mathcal{L} = \frac{1}{{2}} m_1 \dot{u}_1^2+ \frac{1}{{2}} m_2 \dot{u}_2^2-\frac{1}{{2}} \kappa \left( u_2 - u_1 \right)^2,\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

and the pair of Euler-Lagrange equations

\begin{aligned}0 = \frac{d}{dt} \frac{\partial {\mathcal{L}}}{\partial {\dot{u}_i}}- \frac{\partial {\mathcal{L}}}{\partial {u_i}}.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

The equations of motion are

\begin{aligned}\begin{aligned}0 &= m_1 \dot{d}{u}_1 + \kappa \left( u_1 - u_2 \right) \\ 0 &= m_2 \dot{d}{u}_2 + \kappa \left( u_2 - u_1 \right) \end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Let

\begin{aligned}\begin{aligned}u_1(t) &= \int_{-\infty}^\infty \tilde{u}_1(\omega) e^{-i \omega t} d\omega \\ u_2(t) &= \int_{-\infty}^\infty \tilde{u}_2(\omega) e^{-i \omega t} d\omega.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Insertion of these transform pairs into our equations of motion produces a pair of simultaneous integral equations to solve

\begin{aligned}\begin{aligned}0 &= \int_{-\infty}^\infty \left( \left( -m_1 \omega^2 + \kappa \right) \tilde{u}_1(\omega) - \kappa \tilde{u}_2(\omega) \right)e^{-i \omega t} d\omega \\ 0 &= \int_{-\infty}^\infty \left( \left( -m_2 \omega^2 + \kappa \right) \tilde{u}_2(\omega) - \kappa \tilde{u}_1(\omega) \right)e^{-i \omega t} d\omega.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

As with the single spring case, we can decouple these equations with an inverse transformation operation $\int e^{i \omega' t}/2\pi$, which gives us (after dropping primes)

\begin{aligned}0 = \begin{bmatrix}\left( -m_1 \omega^2 + \kappa \right) &- \kappa \\ - \kappa &\left( -m_2 \omega^2 + \kappa \right) \end{bmatrix}\begin{bmatrix}\tilde{u}_1(\omega) \\ \tilde{u}_2(\omega)\end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Taking determinants gives us the constraint on the frequency

\begin{aligned}0 = \left( -m_1 \omega^2 + \kappa \right) \left( -m_2 \omega^2 + \kappa \right) - \kappa^2=m_1 m_2 \omega^4-(m_1 + m_2) \omega^2 =\omega^2 \left( m_1 m_2 \omega^2 -\kappa (m_1 + m_2) \right).\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Introducing a reduced mass

\begin{aligned}\frac{1}{{\mu}} = \frac{1}{m_1}+\frac{1}{m_2},\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

the pair of solutions are

\begin{aligned}\begin{aligned}\omega^2 &= 0 \\ \omega^2 &=\frac{\kappa}{\mu}\equiv \omega_\circ^2.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

As with the single mass oscillator, we require the functions $\tilde{u}{\omega}$ to also be expressed as delta functions. The frequency constraint and that delta function requirement together can be expressed, for $j \in \{0, 1\}$ as

\begin{aligned}\tilde{u}_j(\omega) = A_{j+} \delta( \omega - \omega_\circ )+ A_{j0} \delta( \omega )+ A_{j-} \delta( \omega + \omega_\circ ).\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

With a transformation back to time domain, we have functions of the form

\begin{aligned}\begin{aligned}u_1(t) &= A_{j+} e^{ -i \omega_\circ t }+ A_{j0} + A_{j-} e^{ i \omega_\circ t } \\ u_2(t) &= B_{j+} e^{ -i \omega_\circ t }+ B_{j0} + B_{j-} e^{ i \omega_\circ t }.\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Back insertion of these into the equations of motion, we have

\begin{aligned}\begin{aligned}0 &=-m_1 \omega_\circ^2\left( A_{j+} e^{ -i \omega_\circ t } + A_{j-} e^{ i \omega_\circ t } \right)+ \kappa\left( \left( A_{j+} - B_{j+} \right) e^{ -i \omega_\circ t } + \left( A_{j-} - B_{j-} \right) e^{ i \omega_\circ t } + A_{j0} - B_{j0} \right) \\ 0 &=-m_2 \omega_\circ^2\left( B_{j+} e^{ -i \omega_\circ t } + B_{j-} e^{ i \omega_\circ t } \right)+ \kappa\left( \left( B_{j+} - A_{j+} \right) e^{ -i \omega_\circ t } + \left( B_{j-} - A_{j-} \right) e^{ i \omega_\circ t } + B_{j0} - A_{j0} \right)\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Equality requires identity for all powers of $e^{i \omega_\circ t}$, or

\begin{aligned}\begin{aligned}0 &= B_{j0} - A_{j0} \\ 0 &= -m_1 \omega_\circ^2 A_{j+} + \kappa \left( A_{j+} - B_{j+} \right) \\ 0 &= -m_1 \omega_\circ^2 A_{j-} + \kappa \left( A_{j-} - B_{j-} \right) \\ 0 &= -m_2 \omega_\circ^2 B_{j+} + \kappa \left( B_{j+} - A_{j+} \right) \\ 0 &= -m_2 \omega_\circ^2 B_{j-} + \kappa \left( B_{j-} - A_{j-} \right),\end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

or $B_{j0} = A_{j0}$ and

\begin{aligned}0 =\begin{bmatrix}\kappa -m_1 \omega_\circ^2 & 0 & - \kappa & 0 \\ 0 & \kappa -m_1 \omega_\circ^2 & 0 & - \kappa \\ -\kappa & 0 & \kappa -m_2 \omega_\circ^2 & 0 \\ 0 & -\kappa & 0 & \kappa -m_2 \omega_\circ^2 \end{bmatrix}\begin{bmatrix}A_{j+} \\ A_{j-} \\ B_{j+} \\ B_{j-} \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Observe that

\begin{aligned}\kappa - m_1 \omega_\circ^2=\kappa - m_1 \kappa \left( \frac{1}{{m_1}} + \frac{1}{{m_2}} \right)=\kappa \left( 1 - 1 - \frac{m_1}{m_2} \right)= -\kappa \frac{m_1}{m_2},\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

(with a similar alternate result). We can rewrite eq. 1.0.35.35 as

\begin{aligned}0 =-\kappa\begin{bmatrix}m_1/m_2 & 0 & 1 & 0 \\ 0 & m_1/m_2 & 0 & 1 \\ 1 & 0 & m_2/m_1 & 0 \\ 0 & 1 & 0 & m_2/m_1 \end{bmatrix}\begin{bmatrix}A_{j+} \\ A_{j-} \\ B_{j+} \\ B_{j-} \end{bmatrix}.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

It’s clear that there’s two pairs of linear dependencies here, so the determinant is zero as expected. We can read off the remaining relations. Our undetermined coefficients are given by

\begin{aligned}\begin{aligned}B_{j0} &= A_{j0} \\ m_1 A_{j+} &= -m_2 B_{j+} \\ m_1 A_{j-} &= -m_2 B_{j-} \end{aligned}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

\begin{aligned}\boxed{\begin{aligned}u_1(t) &= a+ m_2 b e^{ -i \omega_\circ t }+ m_2 c e^{ i \omega_\circ t } \\ u_2(t) &= a- m_1 b e^{ -i \omega_\circ t }- m_1 c e^{ i \omega_\circ t }.\end{aligned}}\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

Observe that the constant term is not really of interest, since it represents a constant displacement of both atoms (just a change of coordinates).

Check:

\begin{aligned}u_1(t) - u_2(t) = + (m_1 + m_2)b e^{ -i \omega_\circ t }+ (m_1 + m_2)c e^{ i \omega_\circ t },\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

\begin{aligned}m_1 \dot{d}{u}_1(t) + \kappa (u_1(t) - u_2(t) )=-m_1 m_2 \omega_\circ^2 \left( b e^{ -i \omega_\circ t } + c e^{ i \omega_\circ t } \right)+(m_1 + m_2) \kappa \left( b e^{ -i \omega_\circ t } + c e^{ i \omega_\circ t } \right)=\left( -m_1 m_2 \omega_\circ^2 + (m_1 + m_2) \kappa \right)\left( b e^{ -i \omega_\circ t } + c e^{ i \omega_\circ t } \right)=\left( -m_1 m_2 \kappa \frac{m_1 + m_2}{m_1 m_2} + (m_1 + m_2) \kappa \right)\left( b e^{ -i \omega_\circ t } + c e^{ i \omega_\circ t } \right)= 0.\end{aligned} \hspace{\stretch{1}}(1.0.35.35)

# Reflection

We’ve seen that we can solve any of these constant coefficient systems exactly using matrix methods, however, these will not be practical for large systems unless we have methods to solve for all the non-zero eigenvalues and their corresponding eigenvectors. With the Fourier transform methods we find that our solutions in the frequency domain is of the form

\begin{aligned}\tilde{u}_j(\omega) = \sum a_{kj} \delta( \omega - \omega_{kj} ),\end{aligned} \hspace{\stretch{1}}(1.63)

or in the time domain

\begin{aligned}u(t) = \sum a_{kj} e^{ - i \omega_{kj} t }.\end{aligned} \hspace{\stretch{1}}(1.64)

We assumed exactly this form of solution in class. The trial solution that we used in class factored out a phase shift from $a_{kj}$ of the form of $e^{ i q x_n }$, but that doesn’t change the underling form of that assumed solution. We have however found a good justification for the trial solution we utilized.

## A final pre-exam update of my notes compilation for ‘PHY452H1S Basic Statistical Mechanics’, Taught by Prof. Arun Paramekanti

Posted by peeterjoot on April 22, 2013

Here’s my third update of my notes compilation for this course, including all of the following:

April 21, 2013 Fermi function expansion for thermodynamic quantities

April 20, 2013 Relativistic Fermi Gas

April 10, 2013 Non integral binomial coefficient

April 10, 2013 energy distribution around mean energy

April 09, 2013 Velocity volume element to momentum volume element

April 04, 2013 Phonon modes

April 03, 2013 BEC and phonons

April 03, 2013 Max entropy, fugacity, and Fermi gas

April 02, 2013 Bosons

April 02, 2013 Relativisitic density of states

March 28, 2013 Bosons

plus everything detailed in the description of my previous update and before.

## PHY452H1S Basic Statistical Mechanics. Lecture 21: Phonon modes. Taught by Prof. Arun Paramekanti

Posted by peeterjoot on April 4, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

# Disclaimer

Peeter’s lecture notes from class. May not be entirely coherent.

These are notes for the last class, which included a lot of discussion not captured by this short set of notes (as well as slides which were not transcribed).

# Phonon modes

If we model a solid as a set of interconnected springs, as in fig. 1.1, then the potentials are of the form

Fig 1.1: Solid oscillator model

\begin{aligned}V = \frac{1}{{2}} C \sum_n \left( {u_n - u_{n+1}} \right)^2,\end{aligned} \hspace{\stretch{1}}(1.2.1)

with kinetic energies

\begin{aligned}K = \sum_n \frac{p_n^2}{2m}.\end{aligned} \hspace{\stretch{1}}(1.2.2)

It’s possible to introduce generalized forces

\begin{aligned}F = -\frac{\partial {V}}{\partial {u_n}}\end{aligned} \hspace{\stretch{1}}(1.2.3)

Can differentiate

\begin{aligned}m \frac{d^2 u_n}{dt^2} = - C \left( { u_n - u_{n+1}} \right)- C \left( { u_n - u_{n-1}} \right)\end{aligned} \hspace{\stretch{1}}(1.2.4)

Assuming a Fourier representation

\begin{aligned}u_n = \sum_k \tilde{u}(k) e^{i k n a},\end{aligned} \hspace{\stretch{1}}(1.2.5)

we find

\begin{aligned}m \frac{d^2 \tilde{u}(k)}{dt^2} = - 2 C \left( { 1 - \cos k a} \right)\tilde{u}(k)\end{aligned} \hspace{\stretch{1}}(1.2.6)

This looks like a harmonic oscillator with

\begin{aligned}\omega(k) = \sqrt{ \frac{2 C}{m} ( 1 - \cos k a)}.\end{aligned} \hspace{\stretch{1}}(1.2.7)

This is plotted in fig. 1.2. In particular note that for for $k a \ll 1$ we can use a linear approximation

Fig 1.2: Angular frequency of solid oscillator model

\begin{aligned}\omega(k) \approx \sqrt{ \frac{C}{m} a^2 } \left\lvert {k} \right\rvert.\end{aligned} \hspace{\stretch{1}}(1.2.8)

Experimentally, looking at specific for a complex atomic structure like Gold, we find for example good fit for a model such as

\begin{aligned}C \sim \underbrace{A T}_{\text{Contribution due to electrons.}}+ \underbrace{B T^3}_{\text{Contribution due to phonon like modes where there are linear energy momenum relations.}}.\end{aligned} \hspace{\stretch{1}}(1.2.9)

## An updated compilation of notes, for ‘PHY452H1S Basic Statistical Mechanics’, Taught by Prof. Arun Paramekanti

Posted by peeterjoot on March 3, 2013

That compilation now all of the following too (no further updates will be made to any of these) :

February 28, 2013 Rotation of diatomic molecules

February 28, 2013 Helmholtz free energy

February 26, 2013 Statistical and thermodynamic connection

February 24, 2013 Ideal gas

February 16, 2013 One dimensional well problem from Pathria chapter II

February 15, 2013 1D pendulum problem in phase space

February 14, 2013 Continuing review of thermodynamics

February 13, 2013 Lightning review of thermodynamics

February 11, 2013 Cartesian to spherical change of variables in 3d phase space

February 10, 2013 n SHO particle phase space volume

February 10, 2013 Change of variables in 2d phase space

February 10, 2013 Some problems from Kittel chapter 3

February 07, 2013 Midterm review, thermodynamics

February 06, 2013 Limit of unfair coin distribution, the hard way

February 05, 2013 Ideal gas and SHO phase space volume calculations

February 03, 2013 One dimensional random walk

February 02, 2013 1D SHO phase space

February 02, 2013 Application of the central limit theorem to a product of random vars

January 31, 2013 Liouville’s theorem questions on density and current

January 30, 2013 State counting

## Some problems from Kittel chapter 3

Posted by peeterjoot on February 10, 2013

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

## Question: Classical gas partition function

[1] expresses the classical gas partition function (3.77) as

\begin{aligned}Z_1 \propto \int \exp\left( - \frac{p_x^2 + p_y^2 + p_z^2 }{2 M \tau}\right) dp_x dp_y dp_z\end{aligned} \hspace{\stretch{1}}(1.0.1)

Show that this leads to the expected $3 \tau/2$ result for the thermal average energy.

Let’s use the adjustment technique from the text for the $N$ partition case and write

\begin{aligned}Z_N = \frac{1}{{N!}} Z_1^N,\end{aligned} \hspace{\stretch{1}}(1.0.2)

with $Z_1$ as above. This gives us

\begin{aligned}U &= \tau^2 \frac{\partial {}}{\partial {\tau}} \ln Z_N \\ &= \tau^2 \frac{\partial {}}{\partial {\tau}} \left(N \ln Z_1 - \ln N!\right) \\ &= N \tau^2 \frac{\partial {\ln Z_1 }}{\partial {\tau}} \\ &= N \tau^2 \frac{\partial {}}{\partial {\tau}} \sum_{k = 1}^{3} \ln\int \exp\left( - \frac{p_k^2 }{2 M \tau}\right) dp_k \\ &= N \tau^2\sum_{k = 1}^{3}\frac{\frac{\partial {}}{\partial {\tau}} \int \exp\left( - \frac{p_k^2 }{2 M \tau} \right) dp_k}{\int \exp\left( - \frac{p_k^2 }{2 M \tau} \right) dp_k} \\ &= N \tau^2\sum_{k = 1}^{3}\frac{\frac{\partial {}}{\partial {\tau}} \sqrt{ 2 \pi M \tau }}{\sqrt{ 2 \pi M \tau}} \\ &= 3 N \tau^2\frac{\frac{1}{{2}} \tau^{-1/2}}{\sqrt{ \tau}} \\ &= \frac{3}{2} N \tau \\ &= \frac{3}{2} N k_{\mathrm{B}} T\end{aligned} \hspace{\stretch{1}}(1.0.3)

## Question: Two state system

[1] problem 3.1.

Find an expression for the free energy as a function of $\tau$ of a system with two states, one at energy $0$ and one at energy $\epsilon$. From the free energy, find expressions for the energy and entropy of the system.

Our partition function is

\begin{aligned}Z = 1 + e^{-\epsilon /\tau}\end{aligned} \hspace{\stretch{1}}(1.0.4)

The free energy is just

\begin{aligned}F = -\tau \ln Z = -\tau \ln (1 + e^{-\epsilon/\tau})\end{aligned} \hspace{\stretch{1}}(1.0.5)

The entropy follows immediately

\begin{aligned}\sigma \\ &= -\frac{\partial {F}}{\partial {\tau}} \\ &= \frac{\partial {}}{\partial {\tau}}\left( \tau \ln \left( 1 + e^{-\epsilon/\tau} \right) \right) \\ &= \ln \left( 1 + e^{-\epsilon/\tau} \right)-\tau \epsilon \frac{-1}{\tau^2} \frac{1}{{1 + e^{-\epsilon/\tau}}} \\ &= \ln \left( 1 + e^{-\epsilon/\tau} \right)+\frac{\epsilon}{\tau} \frac{e^{-\epsilon/\tau}}{1 + e^{-\epsilon/\tau}}\end{aligned} \hspace{\stretch{1}}(1.0.6)

The energy is

\begin{aligned}U \\ &= F + \tau \sigma \\ &= -\tau \ln (1 + e^{-\epsilon/\tau}) + \tau \sigma \\ &= \tau\left( \not{{\ln \left( 1 + e^{-\epsilon/\tau} \right)}} + \frac{\epsilon}{\tau} \frac{e^{-\epsilon/\tau}}{1 + e^{-\epsilon/\tau}} -\not{{\ln (1 + e^{-\epsilon/\tau}) }} \right)\end{aligned} \hspace{\stretch{1}}(1.0.7)

This is

\begin{aligned}U=\frac{\epsilon e^{-\epsilon/\tau}}{1 + e^{-\epsilon/\tau}}=\frac{\epsilon}{1 + e^{\epsilon/\tau}}.\end{aligned} \hspace{\stretch{1}}(1.0.8)

These are all plotted in (Fig 1).

Fig1: Plots for two state system

\imageFigure{kittelCh3Problem1PlotsFig1}{Plots for two state system}{fig:kittelCh3Problem1Plots:kittelCh3Problem1PlotsFig1}{0.2}

## Question: Magnetic susceptibility

[1] problem 3.2.

Use the partition function to find an exact expression for the magnetization $M$ and the susceptibility $\chi = dM/dB$ as a function of temperature and magnetic field for the model system of magnetic moments in a magnetic field. The result for the magnetization, found by other means, was $M = n m \tanh( m B/\tau)$, where $n$ is the particle concentration. Find the free energy and express the result as a function only of $\tau$ and the parameter $x = M/nm$. Show that the susceptibility is $\chi = n m^2/\tau$ in the limit $m B \ll \tau$.

Our partition function for a unit volume containing $n$ spins is

\begin{aligned}Z=\frac{\left( e^{-m B/\tau} +e^{m B/\tau} \right)^n}{n!}=2 \frac{\left( \cosh\left( m B/\tau \right) \right)^n}{n!},\end{aligned} \hspace{\stretch{1}}(1.0.9)

so that the Free energy is

\begin{aligned}F = -\tau\left( \ln 2 - \ln n! + n \ln \cosh\left( m B/\tau \right) \right).\end{aligned} \hspace{\stretch{1}}(1.0.15)

The energy, magnetization and magnetic field were interrelated by

\begin{aligned}- M B &= U \\ &= \tau^2 \frac{\partial {}}{\partial {\tau}}\left( -\frac{F}{\tau} \right) \\ &= \tau^2 n\frac{\partial {}}{\partial {\tau}}\ln \cosh\left( m B/\tau \right) \\ &= \tau^2 n \frac{ -m B/\tau^2\sinh\left( m B/\tau \right)}{\cosh\left( m B/\tau \right)} \\ &= - m B n \tanh \left( m B/\tau \right).\end{aligned} \hspace{\stretch{1}}(1.0.11)

This gives us

\begin{aligned}M = m n \tanh \left( m B/\tau \right),\end{aligned} \hspace{\stretch{1}}(1.0.12)

so that

\begin{aligned}\chi = \frac{dM}{dB}= \frac{m^2 n}{\tau \cosh^2 \left( m B/\tau \right)}.\end{aligned} \hspace{\stretch{1}}(1.0.13)

For $m B/\tau \ll 1$, the cosh term goes to unity, so we have

\begin{aligned}\chi \approx= \frac{m^2 n}{\tau},\end{aligned} \hspace{\stretch{1}}(1.0.14)

as desired.

With $x = M/nm$, or $m = M/nx$, the free energy is

\begin{aligned}F =-\tau\left( \ln 2/n! + n \ln \cosh\left( \frac{M B}{n x \tau} \right) \right)\end{aligned} \hspace{\stretch{1}}(1.0.15)

That last expression isn’t particularly illuminating. What was the point of that substitution?

## Question: Free energy of a harmonic oscillator

[1] problem 3.3.

A one dimensional harmonic oscillator has an infinite series of equally spaced energy states, with $\epsilon_s = s \hbar \omega$, where $s$ is a positive integer or zero, and $\omega$ is the classical frequency of the oscillator. We have chosen the zero of energy at the state $s = 0$. Show that for a harmonic oscillator the free energy is

\begin{aligned}F = \tau \ln\left( 1 - e^{-\hbar \omega/\tau} \right).\end{aligned} \hspace{\stretch{1}}(1.0.16)

Note that at high temperatures such that $\tau \gg \hbar \omega$ we may expand the argument of the logarithm to obtain $F \approx \tau \ln (\hbar \omega/\tau)$. From 1.0.16 show that the entropy is

\begin{aligned}\sigma = \frac{\hbar\omega/\tau}{e^{\hbar \omega/\tau} - 1} -\ln\left( 1 - e^{-\hbar \omega/\tau} \right)\end{aligned} \hspace{\stretch{1}}(1.0.17)

I found it curious that this problem dropped the factor of $\hbar\omega/2$ from the energy. Including it we have

\begin{aligned}\epsilon_s = \left( s + \frac{1}{{2}} \right) \hbar \omega,\end{aligned} \hspace{\stretch{1}}(1.0.18)

So that the partition function is

\begin{aligned}Z= \sum_{s = 0}^\infty e^{-\left( s + \frac{1}{{2}} \right) \hbar \omega/\tau }=e^{-\hbar \omega/2\tau}\sum_{s = 0}^\infty e^{-s \hbar \omega/\tau}.\end{aligned} \hspace{\stretch{1}}(1.0.19)

The free energy is

\begin{aligned}F &= -\tau \ln Z \\ &= -\tau\left( -\frac{\hbar \omega}{2\tau} + \ln \left( \sum_{s = 0}^\infty e^{-s \hbar \omega/\tau} \right) \right) \\ &= \frac{\hbar \omega}{2} +\ln\left( \sum_{s = 0}^\infty e^{-s \hbar \omega/\tau} \right)\end{aligned} \hspace{\stretch{1}}(1.0.20)

We see that the contribution of the $\hbar \omega/2$ in the energy of each state just adds a constant factor to the free energy. This will drop out when we compute the entropy. Dropping that factor now that we know why it doesn’t contribute, we can complete the summation, so have, by inspection

\begin{aligned}F = -\tau \ln Z=\tau \ln\left( 1 - e^{-\hbar \omega/\tau} \right).\end{aligned} \hspace{\stretch{1}}(1.0.21)

Taking derivatives for the entropy we have

\begin{aligned}\sigma &= -\frac{\partial {F}}{\partial {\tau}} \\ &= -\ln\left( 1 - e^{-\hbar \omega/\tau} \right)+\tau\frac{\hbar \omega}{\tau^2} \frac{e^{-\hbar \omega/\tau}}{1 - e^{-\hbar \omega/\tau}} \\ &= -\ln\left( 1 - e^{-\hbar \omega/\tau} \right)+\frac{\frac{\hbar \omega}{\tau}}{e^{\hbar \omega/\tau} - 1}\end{aligned} \hspace{\stretch{1}}(1.0.22)

## Question: Energy fluctuation

[1] problem 3.4.

Consider a system of fixed volume in thermal contact with a reservoir. Show that the mean square fluctuation in the energy of the system is

\begin{aligned}\left\langle{{ (\epsilon - \left\langle{\epsilon}\right\rangle)^2 }}\right\rangle = \tau^2\left( \frac{\partial {U}}{\partial {\tau}} \right)_V\end{aligned} \hspace{\stretch{1}}(1.0.23)

Here $U$ is the conventional symbol for $\left\langle{{\epsilon}}\right\rangle$. Hint: Use the partition function $Z$ to relate ${\partial {U}}/{\partial {t}}$ to the mean square fluctuation. Also, multiply out the term $(\cdots)^2$.

With a probability of finding the system in state $s$ of

\begin{aligned}P_s = \frac{e^{-\epsilon_s/\tau}}{Z}\end{aligned} \hspace{\stretch{1}}(1.0.24)

the average energy is

\begin{aligned}U &= \left\langle{{\epsilon}}\right\rangle \\ &= \sum_s P_s \epsilon_s \\ &= \sum_s \epsilon_s \frac{e^{-\epsilon_s/\tau}}{Z} \\ &= \frac{1}{{Z}} \sum_s \epsilon_s e^{-\epsilon_s/\tau}\end{aligned} \hspace{\stretch{1}}(1.0.25)

So we have

\begin{aligned}\tau^2 \frac{\partial {U}}{\partial {\tau}} \\ &= -\frac{\tau^2}{Z^2} \frac{dZ}{d\tau}\sum_s \epsilon_s e^{-\epsilon_s/\tau}+ \frac{\tau^2}{Z}\sum_s \frac{\epsilon_s^2}{\tau^2} e^{-\epsilon_s/\tau} \\ &= -\frac{\tau^2}{Z^2} \frac{dZ}{d\tau}\sum_s \epsilon_s e^{-\epsilon_s/\tau}+ \frac{1}{{Z}}\sum_s \epsilon_s^2 e^{-\epsilon_s/\tau}.\end{aligned} \hspace{\stretch{1}}(1.0.26)

But

\begin{aligned}\frac{dZ}{d\tau}=\frac{d}{d\tau} \sum_s e^{-\epsilon_s/\tau}=\sum_s \frac{\epsilon_s}{\tau^2} e^{-\epsilon_s/\tau},\end{aligned} \hspace{\stretch{1}}(1.0.27)

giving

\begin{aligned}\tau^2 \frac{\partial {U}}{\partial {\tau}} &= \frac{1}{Z^2}\sum_s \epsilon_s e^{-\epsilon_s/\tau} \sum_s \epsilon_s e^{-\epsilon_s/\tau}+ \frac{1}{{Z}}\sum_s \epsilon_s^2 e^{-\epsilon_s/\tau} \\ &= -\left\langle{{ \epsilon }}\right\rangle^2 + \left\langle{{\epsilon^2}}\right\rangle,\end{aligned} \hspace{\stretch{1}}(1.0.28)

which shows 1.0.23 as desired.

# References

[1] C. Kittel and H. Kroemer. Thermal physics. WH Freeman, 1980.

## PHY456H1F: Quantum Mechanics II. Lecture 8 (Taught by Prof J.E. Sipe). Time dependent pertubation (cont.)

Posted by peeterjoot on October 8, 2011

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

# Disclaimer.

Peeter’s lecture notes from class. May not be entirely coherent.

# Time dependent pertubation.

We’d gotten as far as calculating

\begin{aligned}c_m^{(1)}(\infty) = \frac{1}{{i \hbar}} \boldsymbol{\mu}_{ms} \cdot \mathbf{E}(\omega_{ms})\end{aligned} \hspace{\stretch{1}}(2.1)

where

\begin{aligned}\mathbf{E}(t) = \int \frac{d\omega}{2 \pi} \mathbf{E}(\omega) e^{-i \omega t},\end{aligned} \hspace{\stretch{1}}(2.2)

and

\begin{aligned}\omega_{ms} = \frac{E_m - E_s}{\hbar}.\end{aligned} \hspace{\stretch{1}}(2.3)

Graphically, these frequencies are illustrated in figure (\ref{fig:qmTwoL8fig0FrequenciesAbsorbtionAndEmission})

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{qmTwoL8fig0FrequenciesAbsorbtionAndEmission}
\caption{Positive and negative frequencies.}
\end{figure}

The probability for a transition from $m$ to $s$ is therefore

\begin{aligned}\rho_{m \rightarrow s} = {\left\lvert{ c_m^{(1)}(\infty) }\right\rvert}^2= \frac{1}{{\hbar}}^2 {\left\lvert{\boldsymbol{\mu}_{ms} \cdot \mathbf{E}(\omega_{ms})}\right\rvert}^2\end{aligned} \hspace{\stretch{1}}(2.4)

Recall that because the electric field is real we had

\begin{aligned}{\left\lvert{\mathbf{E}(\omega)}\right\rvert}^2 = {\left\lvert{\mathbf{E}(-\omega)}\right\rvert}^2.\end{aligned} \hspace{\stretch{1}}(2.5)

Suppose that we have a wave pulse, where our field magnitude is perhaps of the form

\begin{aligned}E(t) = e^{-t^2/T^2} \cos(\omega_0 t),\end{aligned} \hspace{\stretch{1}}(2.6)

as illustated with $\omega = 10, T = 1$ in figure (\ref{fig:gaussianWavePacket}).

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{gaussianWavePacket}
\caption{Gaussian wave packet}
\end{figure}

We expect this to have a two lobe Fourier spectrum, with the lobes centered at $\omega = \pm 10$, and width proportional to $1/T$.

For reference, as calculated using Mathematica this Fourier transform is

\begin{aligned}E(\omega) = \frac{e^{-\frac{1}{4} T^2 (\omega_0+\omega )^2}}{2 \sqrt{\frac{2}{T^2}}}+\frac{e^{\omega_0 T^2 \omega -\frac{1}{4} T^2 (\omega_0+\omega )^2}}{2 \sqrt{\frac{2}{T^2}}}\end{aligned} \hspace{\stretch{1}}(2.7)

This is illustrated, again for $\omega_0 = 10, and T=1$, in figure (\ref{fig:FTgaussianWavePacket})

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{FTgaussianWavePacket}
\caption{FTgaussianWavePacket}
\end{figure}

where we see the expected Gaussian result, since the Fourier transform of a Gaussian is a Gaussian.

FIXME: not sure what the point of this was?

# Sudden pertubations.

Given our wave equation

\begin{aligned}i \hbar \frac{d{{}}}{dt} {\lvert {\psi(t)} \rangle} = H(t) {\lvert {\psi(t)} \rangle}\end{aligned} \hspace{\stretch{1}}(3.8)

and a sudden pertubation in the Hamiltonian, as illustrated in figure (\ref{fig:suddenStepHamiltonian})

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{suddenStepHamiltonian}
\caption{Sudden step Hamiltonian.}
\end{figure}

Consider $H_0$ and $H_F$ fixed, and decrease $\Delta t \rightarrow 0$. We can formally integrate 3.8

\begin{aligned}\frac{d{{}}}{dt} {\lvert {\psi(t)} \rangle} = \frac{1}{{i \hbar}} H(t) {\lvert {\psi(t)} \rangle}\end{aligned} \hspace{\stretch{1}}(3.9)

For

\begin{aligned}{\lvert {\psi(t)} \rangle} -{\lvert {\psi(t_0)} \rangle} = \frac{1}{{i \hbar}} \int_{t_0}^t H(t') {\lvert {\psi(t')} \rangle} dt'.\end{aligned} \hspace{\stretch{1}}(3.10)

While this is an exact solution, it is also not terribly useful since we don’t know ${\lvert {\psi(t)} \rangle}$. However, we can select the small interval $\Delta t$, and write

\begin{aligned}{\lvert {\psi(\Delta t/2)} \rangle} ={\lvert {\psi(-\Delta t/2)} \rangle}+ \frac{1}{{i \hbar}} \int_{t_0}^t H(t') {\lvert {\psi(t')} \rangle} dt'.\end{aligned} \hspace{\stretch{1}}(3.11)

Note that we could use the integral kernel iteration technique here and substitute ${\lvert {\psi(t')} \rangle} = {\lvert {\psi(-\Delta t/2)} \rangle}$ and then develop this, to generate a power series with $(\Delta t/2)^k$ dependence. However, we note that 3.11 is still an exact relation, and if $\Delta t \rightarrow 0$, with the integration limits narrowing (provided $H(t')$ is well behaved) we are left with just

\begin{aligned}{\lvert {\psi(\Delta t/2)} \rangle} = {\lvert {\psi(-\Delta t/2)} \rangle}\end{aligned} \hspace{\stretch{1}}(3.12)

Or

\begin{aligned}{\lvert {\psi_{\text{after}}} \rangle} = {\lvert {\psi_{\text{before}}} \rangle},\end{aligned} \hspace{\stretch{1}}(3.13)

provided that we change the Hamiltonian fast enough. On the surface there appears to be no consequences, but there are some very serious ones!

## Example: Harmonic oscillator.

Consider our harmonic oscillator Hamiltonian, with

\begin{aligned}H_0 &= \frac{P^2}{2m} + \frac{1}{{2}} m \omega_0^2 X^2 \\ H_F &= \frac{P^2}{2m} + \frac{1}{{2}} m \omega_F^2 X^2\end{aligned} \hspace{\stretch{1}}(3.14)

Here $\omega_0 \rightarrow \omega_F$ continuously, but very quickly. In effect, we have tightened the spring constant. Note that there are cases in linear optics when you can actually do exactly that.

Imagine that ${\lvert {\psi_{\text{before}}} \rangle}$ is in the ground state of the harmonic oscillator as in figure (\ref{fig:suddenHamiltonianPertubationHO})

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{suddenHamiltonianPertubationHO}
\caption{Harmonic oscillator sudden Hamiltonian pertubation.}
\end{figure}

and we suddenly change the Hamilontian with potential $V_0 \rightarrow V_F$ (weakening the “spring”). Professor Sipe gives us a graphical demo of this, by impersonating a constrained wavefunction with his arms, doing weak chicken-flapping of them. Now with the potential weakended, he wiggles and flaps his arms with more freedom and somewhat chaotically. His “wave function” arms are now bouncing around in the new limiting potential (initally over doing it and then bouncing back).

We had in this case the exact relation

\begin{aligned}H_0 {\lvert {\psi_0^{(0)}} \rangle} = \frac{1}{{2}} \hbar \omega_0 {\lvert {\psi_0^{(0)}} \rangle}\end{aligned} \hspace{\stretch{1}}(3.16)

but we also have

\begin{aligned}{\lvert {\psi_{\text{after}}} \rangle} = {\lvert {\psi_{\text{before}}} \rangle} = {\lvert {\psi_0^{(0)}} \rangle}\end{aligned} \hspace{\stretch{1}}(3.17)

and

\begin{aligned}H_F {\lvert {\psi_n^{(f)}} \rangle} = \frac{1}{{2}} \hbar \omega_F \left( n + \frac{1}{{2}} \right) {\lvert {\psi_n^{(f)}} \rangle}\end{aligned} \hspace{\stretch{1}}(3.18)

So

\begin{aligned}{\lvert {\psi_{\text{after}}} \rangle}&={\lvert {\psi_0^{(0)}} \rangle} \\ &=\sum_n {\lvert {\psi_n^{(f)}} \rangle}\underbrace{\left\langle{{\psi_n^{(f)}}} \vert {{\psi_0^{(0)}}}\right\rangle }_{c_n} \\ &=\sum_n c_n {\lvert {\psi_n^{(f)}} \rangle}\end{aligned}

and at later times

\begin{aligned}{\lvert {\psi(t)^{(f)}} \rangle}&={\lvert {\psi_0^{(0)}} \rangle} \\ &=\sum_n c_n e^{i \omega_n^{(f)} t} {\lvert {\psi_n^{(f)}} \rangle},\end{aligned}

whereas

\begin{aligned}{\lvert {\psi(t)^{(o)}} \rangle}&=e^{i \omega_0^{(0)} t} {\lvert {\psi_0^{(0)}} \rangle},\end{aligned}

So, while the wave functions may be exactly the same after such a sudden change in Hamiltonian, the dynamics of the situation change for all future times, since we now have a wavefunction that has a different set of components in the basis for the new Hamiltonian. In particular, the evolution of the wave function is now significantly more complex.

FIXME: plot an example of this.

FIXME: what does Adiabatic mean in this context. The usage in class sounds like it was just “really slow and gradual”, yet this has a definition “Of, relating to, or being a reversible thermodynamic process that occurs without gain or loss of heat and without a change in entropy”.

This is treated in section 17.5.2 of the text [1].

This is the reverse case, and we now vary the Hamiltonian $H(t)$ very slowly.

\begin{aligned}\frac{d{{}}}{dt} {\lvert {\psi(t)} \rangle} = \frac{1}{{i \hbar}} H(t) {\lvert {\psi(t)} \rangle}\end{aligned} \hspace{\stretch{1}}(4.19)

We first consider only non-degenerate states, and at $t = 0$ write

\begin{aligned}H(0) = H_0,\end{aligned} \hspace{\stretch{1}}(4.20)

and

\begin{aligned}H_0 {\lvert {\psi_s^{(0)}} \rangle} = E_s^{(0)} {\lvert {\psi_s^{(0)}} \rangle}\end{aligned} \hspace{\stretch{1}}(4.21)

Imagine that at each time $t$ we can find the “instantaneous” energy eigenstates

\begin{aligned}H(t) {\lvert {\hat{\psi}_s(t)} \rangle} = E_s(t) {\lvert {\hat{\psi}_s(t)} \rangle} \end{aligned} \hspace{\stretch{1}}(4.22)

These states do not satisfy Schr\”{o}dinger’s equation, but are simply solutions to the eigen problem. Our standard strategy in pertubation is based on analysis of

\begin{aligned}{\lvert {\psi(t)} \rangle} = \sum_n c_n(t) e^{- i \omega_n^{(0)} t} {\lvert {\psi_n^{(0)} } \rangle},\end{aligned} \hspace{\stretch{1}}(4.23)

\begin{aligned}{\lvert {\psi(t)} \rangle} = \sum_n b_n(t) {\lvert {\hat{\psi}_n(t)} \rangle},\end{aligned} \hspace{\stretch{1}}(4.24)

we will expand, not using our initial basis, but instead using the instananeous kets. Plugging into Schr\”{o}dinger’s equation we have

\begin{aligned}H(t) {\lvert {\psi(t)} \rangle} &= H(t) \sum_n b_n(t) {\lvert {\hat{\psi}_n(t)} \rangle} \\ &= \sum_n b_n(t) E_n(t) {\lvert {\hat{\psi}_n(t)} \rangle} \end{aligned}

This was complicated before with matrix elements all over the place. Now it is easy, however, the time derivative becomes harder. Doing that we find

\begin{aligned}i \hbar \frac{d{{}}}{dt} {\lvert {\psi(t)} \rangle}&=i \hbar\frac{d{{}}}{dt} \sum_n b_n(t) {\lvert {\hat{\psi}_n(t)} \rangle} \\ &=i \hbar\sum_n \frac{d{{b_n(t)}}}{dt} {\lvert {\hat{\psi}_n(t)} \rangle} + \sum_n b_n(t) \frac{d{{}}}{dt} {\lvert {\hat{\psi}_n(t)} \rangle} \\ &= \sum_n b_n(t) E_n(t) {\lvert {\hat{\psi}_n(t)} \rangle} \end{aligned}

We bra ${\langle {\hat{\psi}_m(t)} \rvert}$ into this

\begin{aligned}i \hbar\sum_n \frac{d{{b_n(t)}}}{dt} \left\langle{{\hat{\psi}_m(t)}} \vert {{\hat{\psi}_n(t)}}\right\rangle+ \sum_n b_n(t) {\langle {\hat{\psi}_m(t)} \rvert}\frac{d{{}}}{dt} {\lvert {\hat{\psi}_n(t)} \rangle} = \sum_n b_n(t) E_n(t) \left\langle{{\hat{\psi}_m(t)}} \vert {{\hat{\psi}_n(t)}}\right\rangle ,\end{aligned} \hspace{\stretch{1}}(4.25)

and find

\begin{aligned}i \hbar\frac{d{{b_m(t)}}}{dt} + \sum_n b_n(t) {\langle {\hat{\psi}_m(t)} \rvert}\frac{d{{}}}{dt} {\lvert {\hat{\psi}_n(t)} \rangle} = b_m(t) E_m(t) \end{aligned} \hspace{\stretch{1}}(4.26)

If the Hamiltonian is changed very very slowly in time, we can imagine that ${\lvert {\hat{\psi}_n(t)} \rangle}'$ is also changing very very slowly, but we are not quite there yet. Let’s first split our sum of bra and ket products

\begin{aligned}\sum_n b_n(t) {\langle {\hat{\psi}_m(t)} \rvert}\frac{d{{}}}{dt} {\lvert {\hat{\psi}_n(t)} \rangle} \end{aligned} \hspace{\stretch{1}}(4.27)

into $n \ne m$ and $n = m$ terms. Looking at just the $n = m$ term

\begin{aligned}{\langle {\hat{\psi}_m(t)} \rvert}\frac{d{{}}}{dt} {\lvert {\hat{\psi}_n(t)} \rangle} \end{aligned} \hspace{\stretch{1}}(4.28)

we note

\begin{aligned}0 &=\frac{d{{}}}{dt} \left\langle{{\hat{\psi}_m(t)}} \vert {{\hat{\psi}_n(t)}}\right\rangle \\ &=\left( \frac{d{{}}}{dt} {\langle {\hat{\psi}_m(t)} \rvert} \right) {\lvert {\hat{\psi}_m(t)} \rangle} \\ + {\langle {\hat{\psi}_m(t)} \rvert} \frac{d{{}}}{dt} {\lvert {\hat{\psi}_m(t)} \rangle} \\ \end{aligned}

Something plus its complex conjugate equals 0

\begin{aligned}a + i b + (a + i b)^{*} = 2 a = 0 \implies a = 0,\end{aligned} \hspace{\stretch{1}}(4.29)

so ${\langle {\hat{\psi}_m(t)} \rvert} \frac{d{{}}}{dt} {\lvert {\hat{\psi}_m(t)} \rangle}$ must be purely imaginary. We write

\begin{aligned}{\langle {\hat{\psi}_m(t)} \rvert} \frac{d{{}}}{dt} {\lvert {\hat{\psi}_m(t)} \rangle} = -i \Gamma_s(t),\end{aligned} \hspace{\stretch{1}}(4.30)

where $\Gamma_s$ is real.

# References

[1] BR Desai. Quantum mechanics with basic field theory. Cambridge University Press, 2009.

## Curious problem using the variational method to find the ground state energy of the Harmonic oscillator (cont.)

Posted by peeterjoot on October 7, 2011

Picking up where this problem was last abandoned.

[Click here for a PDF of this post with nicer formatting (especially if my latex to wordpress script has left FORMULA DOES NOT PARSE errors.)]

# Calculating the Fourier terms

Our wave function, with $\beta=1$ is plotted in figure (\ref{fig:expMinusBetsAbsX})

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{expMinusBetsAbsX}
\caption{Exponential trial function with absolute exponential die off.}
\end{figure}

The zeroth order fitting using the gaussian exponential is found to be

\begin{aligned}\psi_0(x) = \sqrt{2 \beta} \text{erfc}\left(\frac{\beta }{\sqrt{2} \alpha }\right)e^{- \alpha^2 x^2/2 +\beta^2/(2 \alpha^2)} \end{aligned} \hspace{\stretch{1}}(6.16)

With $\alpha = \beta = 1$, this is plotted in figure (\ref{fig:expMinusBetsAbsXfirstOrderFitting}) and can be seen to match fairly well

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{expMinusBetsAbsXfirstOrderFitting}
\caption{First ten orders, fitting harmonic oscillator wavefunctions to this trial function.}
\end{figure}

The higher order terms get small fast, but we can see in figure (\ref{fig:expMinusBetsAbsXtenthOrderFitting}), where a tenth order fitting is depicted that it would take a number of them to get anything close to the sharp peak that we have in our exponential trial function.

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{expMinusBetsAbsXtenthOrderFitting}
\caption{Tenth order harmonic oscillator wavefunction fitting.}
\end{figure}

Note that the all the brakets of even orders in $n$ with the trial function are zero, which is why the tenth order approximation is only a sum of six terms.

Details for this \href{https://github.com/peeterjoot/physicsplay/blob/master/notes/phy456/gaussian\

The question of interest is why if we can approximate the trial function so nicely (except at the origin) even with just a first order approximation (polynomial times gaussian functions where the polynomials are Hankel functions), and we can get an exact value for the lowest energy state using the first order approximation of our trial function, why do we get garbage if we include enough terms that the peak is sharp? It must therefore be important to consider the origin, but how do we give some meaning to the derivative of the absolute value function? The key, supplied when asking in office hours for the course, is to express the absolute value function in terms of Heavyside step functions, for which the derivative can be identified as the delta function.

# Correcting, treating the origin this way.

Here’s how we can express the absolute value function using the Heavyside step

\begin{aligned}{\left\lvert{x}\right\rvert} = x \theta(x) - x \theta(-x),\end{aligned} \hspace{\stretch{1}}(7.17)

where the step function is zero for $x 0$ as plotted in figure (\ref{fig:stepFunction}).

\begin{figure}[htp]
\centering
\includegraphics[totalheight=0.2\textheight]{stepFunction}
\caption{stepFunction}
\end{figure}

Expressed this way, with the identification $\theta'(x) = \delta(x)$, we have for the derivative of the absolute value function

\begin{aligned}{\left\lvert{x}\right\rvert}' &= x' \theta(x) - x' \theta(-x) + x \theta'(x) - x \theta'(-x) \\ &= \theta(x) - \theta(-x) + x \delta(x) + x \delta(-x) \\ &= \theta(x) - \theta(-x) + x \delta(x) + x \delta(x) \\ \end{aligned}

Observe that we have our expected unit derivative for $x > 0$, and $-1$ derivative for $x < 0$. At the origin our $\theta$ contributions vanish, and we are left with

\begin{aligned}{\left.{{{\left\lvert{x}\right\rvert}' }}\right\vert}_{{x=0}}= 2 {\left.{{x \delta(x)}}\right\vert}_{{x=0}} \\ \end{aligned}

We’ve got zero times infinity here, so how do we give meaning to this? As with any delta functional, we’ve got to apply it to a well behaved (square integrable) test function $f(x)$ and integrate. Doing so we have

\begin{aligned}\int_{-\infty}^\infty dx {\left\lvert{x}\right\rvert}' f(x) &= 2 \int_{-\infty}^\infty dx x \delta(x) f(x) \\ &= 2 (0) f(0)\end{aligned}

This equals zero for any well behaved test function $f(x)$. Since the delta function only picks up the contribution at the origin, we can therefore identify ${\left\lvert{x}\right\rvert}'$ as zero at the origin.

Using the same technique, we can express our trial function in terms of steps

\begin{aligned}\psi = e^{-\beta {\left\lvert{x}\right\rvert}} = \theta(x) e^{-\beta x} + \theta(-x) e^{\beta x}.\end{aligned} \hspace{\stretch{1}}(7.18)

This we can now take derivatives of, even at the origin, and find

\begin{aligned}\psi' &= \theta'(x) e^{-\beta x} + \theta'(-x) e^{\beta x} -\beta \theta(x) e^{-\beta x} + \beta \theta(-x) e^{\beta x} \\ &= \delta(x) e^{-\beta x} - \delta(-x) e^{\beta x} -\beta \theta(x) e^{-\beta x} + \beta \theta(-x) e^{\beta x} \\ &= \not{{\delta(x) e^{-\beta x} - \delta(x) e^{\beta x}}} -\beta \theta(x) e^{-\beta x} + \beta \theta(-x) e^{\beta x} \\ &= \beta \left( -\theta(x) e^{-\beta x} + \theta(-x) e^{\beta x} \right)\end{aligned}

Taking second derivatives we find

\begin{aligned}\psi''&= \beta \left( -\theta'(x) e^{-\beta x} + \theta'(-x) e^{\beta x} +\beta \theta(x) e^{-\beta x} + \beta \theta(-x) e^{\beta x} \right) \\ &= \beta \left( -\delta(x) e^{-\beta x} - \delta(-x) e^{\beta x} +\beta \theta(x) e^{-\beta x} + \beta \theta(-x) e^{\beta x} \right) \\ &= \beta^2 \psi - 2 \beta \delta(x)\end{aligned}

Now application of the Hamiltonian operator on our trial function gives us

\begin{aligned}H \psi = -\frac{\hbar^2}{2m} \left( \beta^2 \psi - 2 \delta(x) \right) + \frac{1}{{2}} m \omega^2 x^2 \psi,\end{aligned} \hspace{\stretch{1}}(7.19)

so

\begin{aligned}{\langle {\psi} \rvert} H {\lvert {\psi} \rangle} &= \int_{-\infty}^\infty \left( -\frac{\hbar^2 \beta^2}{2m} + \frac{1}{{2}} m \omega^2 x^2 \right) e^{-2 \beta {\left\lvert{x}\right\rvert} }+ \frac{\hbar^2 \beta}{m}\int_{-\infty}^\infty \delta(x)e^{- \beta {\left\lvert{x}\right\rvert}} \\ &=-\frac{\beta \hbar^2}{2m} + \frac{m \omega^2}{4 \beta^3} + \frac{\hbar^2 \beta}{m} \\ &=\frac{\beta \hbar^2}{2m} + \frac{m \omega^2}{4 \beta^3}.\end{aligned}

Normalized we have

\begin{aligned}E[\beta] = \frac{{\langle {\psi} \rvert} H {\lvert {\psi} \rangle}}{\left\langle{{\psi}} \vert {{\psi}}\right\rangle} = \frac{\beta^2 \hbar^2}{2m} + \frac{m \omega^2}{4 \beta^2}.\end{aligned} \hspace{\stretch{1}}(7.20)

This is looking much more promising. We’ll have the sign alternation that we require to find a positive, non-complex, value for $\beta$ when $E[\beta]$ is minimized. That is

\begin{aligned}0 = \frac{\partial {E}}{\partial {\beta}} = \frac{\beta \hbar^2}{m} - \frac{m \omega^2}{2 \beta^3},\end{aligned} \hspace{\stretch{1}}(7.21)

so the extremum is found at

\begin{aligned}\beta^4 = \frac{m^2 \omega^2}{2 \hbar^2}.\end{aligned} \hspace{\stretch{1}}(7.22)

Plugging this back in we find that our trial function associated with the minimum energy (unnormalized still) is

\begin{aligned}\psi = e^{-\sqrt{\frac{m \omega x^2}{\sqrt{2} \hbar}}},\end{aligned} \hspace{\stretch{1}}(7.23)

and that energy, after substitution, is

\begin{aligned}E[\beta_{\text{min}}] = \frac{\hbar \omega}{2} \sqrt{2}\end{aligned} \hspace{\stretch{1}}(7.24)

We have something that’s $1.4 \times$ the true ground state energy, but is at least a ball park value. However, to get this result, we have to be very careful to treat our point of singularity. A derivative that we’d call undefined in first year calculus, is not only defined, but required, for this treatment to work!

## My submission for PHY356 (Quantum Mechanics I) Problem Set 5.

Posted by peeterjoot on December 7, 2010

[Click here for a PDF of this post with nicer formatting]

The pdf version above has been adjusted after seeing the comments from the grading. [Click here for the PDF for the original submission, as found below, uncorrected.

# Problem.

## Statement

A particle of mass m moves along the x-direction such that $V(X)=\frac{1}{{2}}KX^2$. Is the state

\begin{aligned}u(\xi) = B \xi e^{+\xi^2/2},\end{aligned} \hspace{\stretch{1}}(2.1)

where $\xi$ is given by Eq. (9.60), $B$ is a constant, and time $t=0$, an energy eigenstate of the system? What is probability per unit length for measuring the particle at position $x=0$ at $t=t_0>0$? Explain the physical meaning of the above results.

## Solution

### Is this state an energy eigenstate?

Recall that $\xi = \alpha x$, $\alpha = \sqrt{m\omega/\hbar}$, and $K = m \omega^2$. With this variable substitution Schr\”{o}dinger’s equation for this harmonic oscillator potential takes the form

\begin{aligned}\frac{d^2 u}{d\xi^2} - \xi^2 u = \frac{2 E }{\hbar\omega} u\end{aligned} \hspace{\stretch{1}}(2.2)

While we can blindly substitute a function of the form $\xi e^{\xi^2/2}$ into this to get

\begin{aligned}\frac{1}{{B}} \left(\frac{d^2 u}{d\xi^2} - \xi^2 u\right)&=\frac{d}{d\xi} \left( 1 + \xi^2 \right) e^{\xi^2/2} - \xi^3 e^{\xi^2/2} \\ &=\left( 2 \xi + \xi + \xi^3 \right) e^{\xi^2/2} - \xi^3 e^{\xi^2/2} \\ &=3 \xi e^{\xi^2/2}\end{aligned}

and formally make the identification $E = 3 \omega \hbar/2 = (1 + 1/2) \omega \hbar$, this isn’t a normalizable wavefunction, and has no physical relevance, unless we set $B = 0$.

By changing the problem, this state could be physically relevant. We’d require a potential of the form

\begin{aligned}V(x) =\left\{\begin{array}{l l}f(x) & \quad \mbox{iflatex x < a$} \\ \frac{1}{{2}} K x^2 & \quad \mbox{if$latex a < x b} \\ \end{array}\right.\end{aligned} \hspace{\stretch{1}}(2.3)

For example, $f(x) = V_1, g(x) = V_2$, for constant $V_1, V_2$. For such a potential, within the harmonic well, a general solution of the form

\begin{aligned}u(x,t) = \sum_n H_n(\xi) \Bigl(A_n e^{-\xi^2/2} + B_n e^{\xi^2/2} \Bigr) e^{-i E_n t/\hbar},\end{aligned} \hspace{\stretch{1}}(2.4)

is possible since normalization would not prohibit non-zero $B_n$ values in that situation. For the wave function to be a physically relevant, we require it to be (absolute) square integrable, and must also integrate to unity over the entire interval.

### Probability per unit length at $x=0$.

We cannot answer the question for the probability that the particle is found at the specific $x=0$ position at $t=t_0$ (that probability is zero in a continuous space), but we can answer the question for the probability that a particle is found in an interval surrounding a specific point at this time. By calculating the average of the probability to find the particle in an interval, and dividing by that interval’s length, we arrive at plausible definition of probability per unit length for an interval surrounding $x = x_0$

\begin{aligned}P = \text{Probability per unit length nearlatex x = x_0} =\lim_{\epsilon \rightarrow 0} \frac{1}{{\epsilon}} \int_{x_0 – \epsilon/2}^{x_0 + \epsilon/2} {\left\lvert{ \Psi(x, t_0) }\right\rvert}^2 dx = {\left\lvert{\Psi(x_0, t_0)}\right\rvert}^2\end{aligned} \hspace{\stretch{1}}(2.5)

By this definition, the probability per unit length is just the probability density itself, evaluated at the point of interest.

Physically, for an interval small enough that the probability density is constant in magnitude over that interval, this probability per unit length times the length of this small interval, represents the probability that we will find the particle in that interval.

### Probability per unit length for the non-normalizable state given.

It seems possible, albeit odd, that this question is asking for the probability per unit length for the non-normalizable $E_1$ wavefunction 2.1. Since normalization requires $B=0$, that probability density is simply zero (or undefined, depending on one’s point of view).

### Probability per unit length for some more interesting harmonic oscillator states.

Suppose we form the wavefunction for a superposition of all the normalizable states

\begin{aligned}u(x,t) = \sum_n A_n H_n(\xi) e^{-\xi^2/2} e^{-i E_n t/\hbar}\end{aligned} \hspace{\stretch{1}}(2.6)

Here it is assumed that the $A_n$ coefficients yield unit probability

\begin{aligned}\int {\left\lvert{u(x,0)}\right\rvert}^2 dx = \sum_n {\left\lvert{A_n}\right\rvert}^2 = 1\end{aligned} \hspace{\stretch{1}}(2.7)

For the impure state of 2.6 we have for the probability density

\begin{aligned}{\left\lvert{u}\right\rvert}^2&=\sum_{m,n}A_n A_m^{*} H_n(\xi) H_m(\xi) e^{-\xi^2} e^{-i (E_n - E_m)t_0/\hbar} \\ &=\sum_n{\left\lvert{A_n}\right\rvert}^2 (H_n(\xi))^2 e^{-\xi^2}+\sum_{m \ne n}A_n A_m^{*} H_n(\xi) H_m(\xi) e^{-\xi^2} e^{-i (E_n - E_m)t_0/\hbar} \\ &=\sum_n{\left\lvert{A_n}\right\rvert}^2 (H_n(\xi))^2 e^{-\xi^2}+\sum_{m \ne n}A_n A_m^{*} H_n(\xi) H_m(\xi) e^{-\xi^2} e^{-i (E_n - E_m)t_0/\hbar} \\ &=\sum_n{\left\lvert{A_n}\right\rvert}^2 (H_n(\xi))^2 e^{-\xi^2}+\sum_{m < n}H_n(\xi) H_m(\xi)\left(A_n A_m^{*}e^{-\xi^2} e^{-i (E_n - E_m)t_0/\hbar}+A_m A_n^{*}e^{-\xi^2} e^{-i (E_m - E_n)t_0/\hbar}\right) \\ &=\sum_n{\left\lvert{A_n}\right\rvert}^2 (H_n(\xi))^2 e^{-\xi^2}+2 \sum_{m < n}H_n(\xi) H_m(\xi)e^{-\xi^2}\Re \left(A_n A_m^{*}e^{-i (E_n - E_m)t_0/\hbar}\right) \\ &=\sum_n{\left\lvert{A_n}\right\rvert}^2 (H_n(\xi))^2 e^{-\xi^2} \\ &\quad+2 \sum_{m < n}H_n(\xi) H_m(\xi)e^{-\xi^2}\left(\Re ( A_n A_m^{*} ) \cos( (n - m)\omega t_0)+\Im ( A_n A_m^{*} ) \sin( (n - m)\omega t_0)\right) \\ \end{aligned}

Evaluation at the point $x = 0$, we have

\begin{aligned}{\left\lvert{u(0,t_0)}\right\rvert}^2=\sum_n{\left\lvert{A_n}\right\rvert}^2 (H_n(0))^2 +2 \sum_{m < n} H_n(0) H_m(0) \left( \Re ( A_n A_m^{*} ) \cos( (n - m)\omega t_0) +\Im ( A_n A_m^{*} ) \sin( (n - m)\omega t_0)\right)\end{aligned} \hspace{\stretch{1}}(2.8)

It is interesting that the probability per unit length only has time dependence for a mixed state.

For a pure state and its wavefunction $u(x,t) = N_n H_n(\xi) e^{-\xi^2/2} e^{-i E_n t/\hbar}$ we have just

\begin{aligned}{\left\lvert{u(0,t_0)}\right\rvert}^2=N_n^2 (H_n(0))^2 = \frac{\alpha}{\sqrt{\pi} 2^n n!} H_n(0)^2\end{aligned} \hspace{\stretch{1}}(2.9)

This is zero for odd $n$. For even $n$ is appears that $(H_n(0))^2$ may equal $2^n$ (this is true at least up to n=4). If that’s the case, we have for non-mixed states, with even numbered energy quantum numbers, at $x=0$ a probability per unit length value of ${\left\lvert{u(0,t_0)}\right\rvert}^2 = \frac{\alpha}{\sqrt{\pi} n!}$.